!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! This module provides all the basic facilities for describing and
!! computing the connectivity of a \f$d\f$-dimensional grid immersed
!! in a \f$\mathbb{R}^m\f$, with \f$d\leq m\f$. This information
!! represents the background for the definition of degrees of freedom
!! in physical space, prescribing boundary conditions and so on.
!!
!! \details
!!
!! \section notation Notation
!!
!! In the following we use \f$\{\ldots\}\f$ to denote a set (a
!! collection of elements without any ordering) and \f$[\ldots]\f$ to
!! denote an ordered tuple.
!!
!! \section simplex_definition Definition of the simplex grid
!! 
!! The first concept used in the definition of the grid is the
!! \f$p\f$-simplex. A \f$p\f$-simplex \f$S\f$ is the convex hull of
!! \f$p+1\f$ vertices \f$v_i\f$, each vertex being a point in the
!! Euclidean space \f$\mathbb{R}^m\f$, for \f$p\leq d\f$. Assuming
!! that a (global) collection of vertices is given, \f$[v_i],\,
!! i=1,\ldots,n_v\f$, \f$S\f$ is identified by its \f$p+1\f$ vertices:
!! \f$\{v^S_1,\ldots,v^S_{p+1}\}\f$, taken in any order, or,
!! equivalently, by the corresponding indexes
!! \f$\{i(v^S_1),\ldots,i(v^S_{p+1})\}\f$, with \f$v_{i(v^S)}=v^S\f$.
!! A \f$p\f$-simplex can be oriented by fixing an ordering for its
!! vertices: \f$[v^S_1,\ldots,v^S_{p+1}]\f$. At this point, any
!! positive permutation of the vertices identifies the oriented
!! simplex, while any negative permutation identifies its negative:
!! \f{displaymath}{
!!  \begin{array}{rcrc}
!!     \pi^+(v^S_1,\ldots,v^S_{p+1}) & \mapsto &  S &
!!           {\rm (positive\,\, orientation)} \\\
!!     \pi^-(v^S_1,\ldots,v^S_{p+1}) & \mapsto & -S & 
!!           {\rm (negative\,\, orientation)}
!!  \end{array}
!! \f}
!! where \f$\pi^+\f$ and \f$\pi^-\f$ are positive and negative
!! permutations\latexonly\footnote{This corresponds to
!! \begin{itemize}
!!   \item fixing the orientation $\omega$ on $\mathbb{R}^p$ as $du^1
!!   \wedge ... \wedge du^p$ and 
!!   \item mapping the canonical simplex $\hat{S}$ with vertices
!!   $(0,...,0),(1,...,0),...,(0,...,1)$ on $S$ preserving the order
!!   of the vertices).
!! \end{itemize}}\endlatexonly. 
!! The boundary \f$\partial S\f$ of a \f$p\f$-simplex \f$S\f$ is
!! composed of \f$p+1\f$ \f$(p-1)\f$-simplexes \f$S_b(S)\f$, which can
!! be univocally associated with the vertices of \f$S\f$ by
!! \f{displaymath}{
!!   S_b\, |\, S_b \subset \partial S   \mapsto   v^S\, 
!!         |\, v^S\,\, {\rm vertex\,\, of\,\, }S,
!!          \, v^S\,\, {\rm not\,\, vertex\,\, of\,\,}S_b
!! \f}
!! (essentially, a face is associated with the opposite vertex).
!! 
!! A \f$d\f$-dimensional grid is a collection of three families of
!! simplexes:\latexonly\newline\endlatexonly
!! \f{tabular}{{rl}
!!       $0$-simplexes: & vertices       \\\
!!   $(d-1)$-simplexes: & sides (facets) \\\
!!       $d$-simplexes: & elements (bodies).
!! \f}
!! 
!! A \f$p\f$-simplex \f$S\f$ (typically, \f$S\f$ will be either an
!! element or a side) can be regarded as a \f$p\f$-dimensional
!! manifolds, on which a metric is induced by the scalar product of
!! \f$\mathbb{R}^m\f$. A map for the simplex is readily obtained in
!! terms of the barycentric coordinates
!! \f$\lambda_S^2,\ldots,\lambda_S^{p+1}\f$ by setting
!! \f{displaymath}{
!!   x = \mathscr{F}_S(\lambda_S^2,\ldots,\lambda_S^{p+1}) =
!!   \sum_{i=2}^{p+1}(v^S_i-v^S_1)\lambda_S^i + v^S_1,
!! \f}
!! where the barycentric coordinates vary in the canonical
!! \f$p\f$-simplex
!! \f{displaymath}{
!!  \hat{S} = \left\{ (\lambda_S^2,\ldots,\lambda_S^{p+1})\in
!!  [0\,,1]^p\quad:\quad \sum_{i=2}^{p+1}\lambda_S^i \leq 1\, \right\}.
!! \f}
!! \note It is convenient to define also
!! \f{displaymath}{
!!  \lambda_S^1=1-\sum_{i=2}^{p+1}\lambda_S^i.
!! \f}
!!
!! The contravariant basis induced by this map is
!! \f{displaymath}{
!!   e^S_i = \frac{\partial \mathscr{F}_S}{\partial \lambda_S^{i+1}} = 
!!   v^S_{i+1}-v^S_1, \qquad i=1,\ldots,p.
!! \f}
!! We can now define the metric tensor on \f$S\f$ and the associated
!! volume element by
!! \f{displaymath}{
!!   g_{ij} = e^S_i \cdot e^S_j
!! \f}
!! and
!! \f{displaymath}{
!!   *(1) = \sqrt{det(g_{ij})} \, d\lambda_S^2\wedge\ldots\wedge
!!   d\lambda_S^{p+1}
!! \f}
!! where \f$*\f$ is the Hodge operator. Since the metric is constant
!! on the simplex, and 
!! \f{displaymath}{
!!  vol(\hat{S}) = \frac{1}{p!},
!! \f}
!! we readily have
!! \f{displaymath}{
!!  vol(S) = \frac{1}{p!} \sqrt{det(g_{ij})}.
!! \f}
!! The outward unit normal on \f$\partial S\f$ is also a very
!! important geometrical object. We define it by transforming the
!! outward normal of the reference simplex (notice that the normal
!! transforms, up to a normalization coefficient, as a
!! <em>covariant</em> vector): first define the covariant basis
!! \f$\varepsilon^i_S\f$ by
!! \f{displaymath}{
!!  \varepsilon^i_S \cdot e^S_j = \delta^i_j
!! \f}
!! then let
!! \f{displaymath}{
!!  n^i = \frac{1}{|\tilde{n}^i|} \tilde{n}^i, \qquad {\rm with}\qquad
!!  \tilde{n}^i = \sum_{j=1}^p \hat{n}^i_j \varepsilon^j_S.
!! \f}
!! Notice that the coefficients \f$\hat{n}^i_j\f$ are readily obtained
!! as (see the details in \c mod_master_el)
!! \f{displaymath}{
!!  \hat{n}^1_j = \frac{\sqrt{p}}{p}, \qquad
!!  \hat{n}^i_j = -\delta^{i-1}_{j}.
!! \f}
!! We discuss now some more specific issues
!!
!! \subsection simplex_orientation Simplex orientation
!!
!! Given a simplex \f$S\f$, each permutation of the simplex vertices
!! defines a different map \f$\mathscr{F}_S:\hat{S}\to S\f$, for a
!! total of \f$(p+1)!\f$ possible maps. Fixing the orientation of the
!! simplex halves the number of possible maps. For elements, we can
!! pick any orientation and store it in the mesh structure without
!! further complications. The situation is not as simple for sides,
!! since each element sharing a side will induce an ordering of the
!! side vertices as follows: we orient a generic side \f$\hat{S}_b\f$
!! of \f$\hat{S}\f$ so that the area element \f$\frac{1}{(d-1)!} *
!! (\hat{e}_1^{\hat{S}_b}\wedge\ldots\wedge\hat{e}_{d-1}^{\hat{S}_b})\f$
!! is pointing inwards and then we choose any ordering of the side
!! vertices which is consistent with this orientation, the map
!! \f$\mathscr{F}_S\f$ then transfers this ordering from
!! \f$\hat{S}_b\f$ to \f$S_b\f$. Thus, in general there will be three
!! different orderings for each side: two of them induced by the
!! connected elements and the third one (possibly coinciding with one
!! of the two previous ones) intrinsic to the side. To consistently
!! deal with all of them, each element stores \f$d+1\f$ pointers to
!! the permutations mapping the side ordering induced by the element
!! itself to the intrinsic ordering of each side in the field
!! <tt>e\%pi(:)</tt>, as well as pointers to the inverse permutations
!! from the intrinsic ordering to the element induced one in the field
!! <tt>e\%ip(:)</tt>.
!!
!! \subsection the_0D_case The zero-dimensional case: grid vertices
!!
!! Vertices have no orientation nor volume. For \f$d=1\f$, sides
!! degenerate into \f$0\f$-simplexes; nevertheless they are not
!! regarded as vertices, and the concepts of surface and orientation
!! are suitably extended to maintain their formal properties. See the
!! additional comments in \c t_s and \c t_e for further details.
!!
!! \subsection matrepressimp Matrix representation of the simplex map
!!
!! The map \f$\mathscr{F}_{S}\f$ is represented as
!! \f{displaymath}{
!!   \mathscr{F}_{S}(\xi) = B^S\,\xi + v^S_1, \qquad
!!   B^S = \left[ e^S_1 \, | \ldots | e^S_p\right],
!! \f}
!! i.e. the columns of the matrix \f$B^S\f$ are the vectors \f$e^S\f$.
!! Notice that \f$ \xi = \left[ \lambda^2\,,\ldots\,,\lambda^{p+1}
!! \right]^T \f$. When \f$p=m\f$, the matrix \f$B^S\f$ is square and
!! invertible; when \f$p<m\f$ the matrix \f$B^S\f$ is rectangular and
!! it is possible to define a right inverse \f$B^{S,-1}\f$ such that
!! \f{displaymath}{
!!  B^S\,B^{S,-1} x = x \qquad \forall x\in S.
!! \f}
!! More precisely, since \f$B^S\f$ has full column rank (unless the
!! simplex is degenerate, which we exclude), we can compute the
!! Moore-Penrose pseudoinverse
!! \f{displaymath}{
!!  B^{S,-1}=((B^S)^TB^S)^{-1} (B^S)^T.
!! \f}
!! Since in practice we don't expect large matrices, the computation
!! of \f$B^{S,-1}\f$ doesn't require very sophisticated algorithms.
!! 
!! \section domain_boundaries Domain boundaries
!!
!! Concerning the boundaries of the computational domain, we adopt the
!! following conventions. The boundary sides are always numbered after
!! the internal ones, so that they correspond to the index range \c
!! grid\%ni+1\f$,\ldots,\f$\c grid\%ns, with \c grid\%nb\f$=\f$\c
!! grid\%ns\f$-\f$<tt>grid\%ni</tt>. Each boundary side is tagged with
!! boundary marker \f$m\in\mathbb{N}^+\f$, which shall be specified in
!! the grid input file. Boundary sides are connected to exactly one
!! element, which index is stored in <tt>grid\%s(ib)\%ie(1)</tt>, for
!! \c ib index of a boundary side, while the boundary marker with
!! changed sign is stored in <tt>grid\%s(ib)\%ie(2)</tt>. The boundary
!! information is also repeated in the boundary elements, in the
!! fields <tt>grid\%e(ieb)\%ie(ibl)</tt> and
!! <tt>grid\%e(ieb)\%e(ibl)\%p</tt>, where \c ieb is the index of the
!! boundary element and \c ibl is the local index of the boundary.
!! Summarizing, boundary sides/elements can be detected by any of the
!! following conditions:
!! - boundary sides
!!   - \c ib \f$>\f$ \c grid\%ni,
!!   - <tt>grid\%s(ib)\%ie(2)</tt>\f$<0\f$,
!!   - <tt>.not.associated(grid\%s(ib)\%e(2)\%p)</tt>,
!!   .
!! - boundary elements
!!   - <tt>grid\%e(ieb)\%ie(ibl)</tt>\f$<0\f$,
!!   - <tt>.not.associated(grid\%e(ieb)\%e(ibl)\%p)</tt>.
!!   .
!! .
!! \note In this module the boundary marker is simply stored in \c
!! grid, there is no distinction between Dirichlet or Neumann boundary
!! condition.
!!
!! \section domain_decomposition Domain decomposition
!!
!! Concerning domain decomposition information, if a side belongs to
!! only one element but is not listed as a boundary side in \c e (see
!! \c new_grid), than it is assumed that it is a domain decomposition
!! side. Such sides are considered boundary sides (so that they comply
!! to all the specifications given in section \ref domain_boundaries)
!! and get assigned a boundary marker larger than the largest marker
!! specified in \c e (notice that this requires a global reduction,
!! since \c e is known only for the local subgrid).  If there is at
!! least one domain decomposition side, the relevant information is
!! read from the grid file and returned by \c new_grid.  Essentially,
!! this information is collected in a \c t_ddc_grid object and is of
!! two types:
!! <ul>
!!  <li> each entity needs a global index which is unique in the
!!   <em>global</em> domain
!!  <li> each shared entity, i.e. each node or side which is shared by
!!   two or more domains, requires some translation information
!!   specifying how the same entity is seen by the connected domains
!!   (notice that, by construction, there are no shared elements).
!! </ul>
!! The first kind of information makes it possible to define a global
!! structure for the mesh, leading, for instance, to a global
!! numbering of the degrees of freedom; this is required in the
!! solution of global linear systems or in the visualization of global
!! fields. The second kind of information allows each domain to
!! communicate with its neighbours; this is required, for instance, in
!! the exchange of numerical fluxes in a finite volume context.
!! \note The additional domain decomposition information in the grid
!! file is required only if there is at least one domain decomposition
!! side.
!<----------------------------------------------------------------------
module mod_grid

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_integer,                  &
   mpi_max,                      &
   mpi_status_size,              &
   mpi_comm_size, mpi_comm_rank, &
   mpi_request_null,             &
   mpi_irecv, mpi_isend,         &
   mpi_bcast,                    &
   mpi_wait, mpi_waitall,        &
   mpi_allreduce,                &
   mpi_alltoall, mpi_alltoallv

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave,   &
   read_octave,    &
   read_octave_al, &
   locate_var

 use mod_perms, only: &
   mod_perms_constructor, &
   mod_perms_destructor,  &
   mod_perms_initialized, &
   t_perm, t_dperm, &
   operator(.eq.),  &
   operator(.ne.),  &
   operator(.lt.),  &
   operator(.le.),  &
   operator(.gt.),  &
   operator(.ge.),  &
   operator(*),     &
   operator(**),    &
   perm_table,      &
   idx,             &
   dperm_reduce,    &
   fact

 use mod_linal, only: &
   mod_linal_initialized, &
   qr, invmat

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me, new_me, clear, &
   write_octave

!-----------------------------------------------------------------------
 
 implicit none
 save

!-----------------------------------------------------------------------
 
! Module interface

 public :: &
   mod_grid_constructor, &
   mod_grid_destructor,  &
   mod_grid_initialized, &
   t_v, t_s, t_e,       &
   p_t_v, p_t_s, p_t_e, &
   t_ddc_gv, t_ddc_nv,  &
   t_grid, t_ddc_grid,  &
   new_grid, clear,     &
   new_subgrid,         &
   affmap, iaffmap,     &
   diameter,            &
   el_linear_size,      & 
   el_courant,          &
   locate_point,        &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 ! pointer arrays
 type p_t_v
   type(t_v), pointer :: p => null()
 end type p_t_v
 type p_t_s
   type(t_s), pointer :: p => null()
 end type p_t_s
 type p_t_e
   type(t_e), pointer :: p => null()
 end type p_t_e

 !> vertex
 type t_v
   !> vertex properties
   integer :: m !< space dimension
   integer :: i !< vertex index
   real(wp), allocatable :: x(:) !< coordinates (m)
   !> connected sides
   integer :: ns                    !< # connected sides
   integer, allocatable :: is(:)    !< ind. connected sides
   type(p_t_s), allocatable :: s(:) !< connected sides
   !> connected elements
   integer :: ne                    !< # connected elements
   integer, allocatable :: ie(:)    !< ind. connected elements
   type(p_t_e), allocatable :: e(:) !< connected elements
 end type t_v

 !> sides
 !!
 !! A side respects the following constraints:
 !! <ol>
 !!  <li> the vertex ordering used in \c iv and \c v defines the side
 !!  orientation, all the geometric quantities are consistent with
 !!  this orientation
 !!  <li> boudary sides have <tt>s\%i.gt.grid\%ni</tt> and
 !!  <tt>s\%ie(2).lt.0</tt> and <tt>s\%e(2)=>null()</tt>.
 !! </ol>
 !! \note For the 1D case, some fields need an <em>ad hoc</em>
 !! definition. These are: \c xb and \c x0 coincide with the point
 !! representing the side, \c a is always 1 and \c b is a
 !! \f$m\times0\f$ array. All the functions provided in this module
 !! deal correctly with the 1D case.
 !!
 !! \bug The affine map should be a type bound procedure
 !<
 type t_s
   integer :: d, m                  !< dimensions
   integer :: i                     !< side index
   integer, allocatable :: iv(:)    !< ind. vertices (d)
   type(p_t_v), allocatable :: v(:) !< vertices      (d)
   !> connected elements
   integer :: ie(2)                 !< ind. elements
   type(p_t_e) :: e(2)              !< elements
   integer :: isl(2)                !< local indx. on e
   !> geometry
   real(wp), allocatable :: xb(:)   !< barycenter (m)
   real(wp) :: a                    !< surface
   !> affine map: \f$x = b\,\xi + x_0\f$
   real(wp), allocatable :: b(:,:)  !< (m,d-1)
   real(wp), allocatable :: x0(:)   !< (m)
 end type t_s

 !> element
 !!
 !! \note The field \c det_b is in fact computed as
 !! \f$\displaystyle\frac{vol(E)}{vol(\hat{E})}\f$. When \f$d=m\f$,
 !! \f$B^E\f$ is square and this corresponds to \f$det(B^E)\f$.
 !! \bug the affine map and the outward normal should be type bound
 !! procedures
 type t_e
   integer :: d, m !< dimensions
   integer :: i    !< element index
   integer, allocatable :: iv(:)    !< ind. vertices (d+1)
   type(p_t_v), allocatable :: v(:) !< vertices      (d+1)
   !> connected sides
   integer, allocatable :: is(:)    !< ind. sides    (d+1)
   type(p_t_s), allocatable :: s(:) !< sides         (d+1)
   integer, allocatable :: pi(:)    !< el2side perm. (d+1)
   integer, allocatable :: ip(:)    !< side2el perm. (d+1)
   !> connected elements
   integer, allocatable :: ie(:)    !< ind. neigh. elements (d+1)
   type(p_t_e), allocatable :: e(:) !< neighboring elements (d+1)
   integer, allocatable :: iel(:)   !< local indx. on neigs (d+1)
   !> geometry
   real(wp), allocatable :: xb(:)   !< barycenter (m)
   real(wp) :: vol                  !< volume
   !> affine map: \f$x = b\,\xi + x_0\f$, inverse \f$\xi = b_i(x-x_0)\f$
   real(wp), allocatable :: b(:,:)  ! (m,d)
   real(wp) :: det_b                ! det(b)
   real(wp), allocatable :: bi(:,:) ! (d,m)
   real(wp), allocatable :: x0(:)   ! (m)
   real(wp), allocatable :: n(:,:)  ! normals (m,d+1)
   !> Additional data
   !!
   !! This field can be used to store any data which are attached to
   !! the elements. These data are simply copied from the grid file,
   !! they don't have to have any specific meaning here. A field \c
   !! edata_legend is included in \c t_grid to provide a description
   !! of these values.
   real(wp), allocatable :: edata(:)
 end type t_e

 !> grid
 !> \note all the variables of type t_grid should have the TARGET
 !> attribute, so that each element can point to the permutation table
 !! and also all the internal pointers can be set.
 !> \bug the fields d, m should become LEN parameters
 type t_grid
   integer :: d, m !< dimensions
   integer :: ne !< # elements
   integer :: nv !< # vertices
   integer :: ns !< ns = ni+nb
   integer :: nb !< # boundary sides
   integer :: ni !< # internal sides
   type(t_me) :: me !< master element
   type(t_v), allocatable :: v(:) !< vertices
   type(t_s), allocatable :: s(:) !< sides
   type(t_e), allocatable :: e(:) !< elements
   !> geometry
   real(wp) :: vol !< total volume
   !> Additional data (see <tt>t_e\%edata</tt>)
   character(len=100), allocatable :: edata_legend(:)
 end type t_grid

 !> domain decomposition global vertex
 type t_ddc_gv
   integer :: gi !< global vertex index
   !> neighbour vertex index, larger than zero if the vertex belongs
   !! to more than one domain, in which case it gives the index of the
   !! \c t_ddc_nv corresponding object
   integer :: ni = 0
 end type t_ddc_gv
 !> domain decomposition neighbour vertex
 type t_ddc_nv
   integer :: i  !< local vertex index
   integer :: nd !< number of neighboring domains sharing the vertex
   integer, allocatable :: id(:) !< neighboring domains (nd)
   integer, allocatable :: in(:) !< neighbour local indexes (nd)
   !> interface index
   !!
   !! For each pair of subdomains, the shared vertexes receive a
   !! unique interface index, which is then useful for the pairwise
   !! communications.
   !!
   !! \todo Implement the interface index
   !<
   !integer, allocatable :: i_itf(:)
 end type t_ddc_nv

 !> domain decomposition global side
 type t_ddc_gs
   integer :: gi !< global side index
   !> neighbour side index, larger than zero if the side belongs to
   !! more than one domain, in which case it gives the index of the \c
   !! t_ddc_ns corresponding object
   !<
   integer :: ni = 0
 end type t_ddc_gs
 !> domain decomposition neighbour side
 type t_ddc_ns
   integer :: i  !< local side index
   integer :: id !< neighboring domain
   integer :: in !< neighbour local index
   !> side to side permutation
   !!
   !! This is the permutation required to map the vertices from this
   !! domain (side ordering) to the neighboring domain (side ordering)
   !<
   integer :: p_s2s
   !> interface index
   !!
   !! For each pair of subdomains, the shared sides receive a
   !! unique interface index, which is then useful for the pairwise
   !! communications.
   !!
   !! \todo Implement the interface index
   !<
   !integer :: i_itf
 end type t_ddc_ns

 !> domain decomposition global element
 type t_ddc_ge
   integer :: gi !< global element index
 end type t_ddc_ge

 !> domain decomposition grid
 !!
 !! \note The indexes in \c nnv_id and \c nns_id start from zero, so
 !! that these arrays can be indexed directly with the domain id-s.
 !! \note A neighbour vertex can be counted for more than one
 !! neighboring domain in \c nnv_id, so that in general
 !! <tt>sum(nnv_id).ge.nnv</tt>. This is not the case for neighbour
 !! sides.
 !<
 type t_ddc_grid
   integer :: nd !< number of subdomains
   integer :: id !< subdomain index (starts from 0)
   integer :: ngv !< total number of vertexes
   integer :: ngs !< total number of sides
   integer :: nge !< total number of elements
   integer :: nnv !< \# neighbour vertices
   integer :: nns !< \# neighbour sides
   !> number of neighbour vertices for each domain (0:nd-1)
   integer, allocatable :: nnv_id(:)
   !> number of neighbour sides for each domain (0:nd-1)
   integer, allocatable :: nns_id(:)
   integer :: ddc_marker !< marker for the domain decomposition sides
   type(t_ddc_gv), allocatable :: gv(:) !< global vertices (nv)
   type(t_ddc_gs), allocatable :: gs(:) !< global sides (ns)
   type(t_ddc_ge), allocatable :: ge(:) !< global elements (ne)
   type(t_ddc_nv), allocatable :: nv(:) !< neighbour vertices (nnv)
   type(t_ddc_ns), allocatable :: ns(:) !< neighbour sides (nns)
 end type t_ddc_grid

 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_grid_initialized = .false.

 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_grid'

 interface new_grid
   module procedure new_grid, new_grid_file
 end interface

 interface clear
   module procedure clear_grid, clear_ddc_grid
 end interface

 interface affmap
   module procedure affmap, affmap_s
 end interface

 interface iaffmap
   module procedure iaffmap
 end interface

 interface diameter
   module procedure diameter_e, diameter_s
 end interface

 interface write_octave
   module procedure write_grid_struct,     write_vert_struct,   &
                    write_side_struct,     write_elem_struct,   &
                    write_ddc_grid_struct, write_ddc_vert_struct
 end interface write_octave

 interface read_octave_al
   module procedure read_al_ddc_v_struct
 end interface read_octave_al

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_grid_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_mpi_utils_initialized.eqv..false.)  .or. &
       (mod_octave_io_initialized.eqv..false.)  .or. &
       (mod_perms_initialized.eqv..false.)      .or. &
       (mod_master_el_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_grid_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_grid_initialized = .true.
 end subroutine mod_grid_constructor

!-----------------------------------------------------------------------

 subroutine mod_grid_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_grid_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_grid_initialized = .false.
 end subroutine mod_grid_destructor

!-----------------------------------------------------------------------
 
 !> Build a \c t_grid object reading data from file.
 !! 
 !! The input file shall contain four variables:
 !! <ul>
 !!  <li> d: dimension of the grid (scalar)
 !!  <li> p: coordinates of the vertices (m,nv)
 !!  <li> e: side information (d+3,nb) (see later)
 !!  <li> t: element connectivity (d+1,ne) (additional rows ignored)
 !! </ul>
 !! For additional details see \c new_grid. If the grid is obtained by
 !! domain decomposition, the related information must be also
 !! provided in the grid file as described in \c new_grid.
 !!
 !! If the grid file contains a character array \c edata_legend, this
 !! is copied in <tt>grid\%edata_legend</tt> and the corresponding
 !! element data are read from an array \c edata, which should have
 !! one column for each element. If \c edata_legend is not present, the
 !! corresponding edata fields are left unallocated.
 subroutine new_grid_file( grid, grid_file_name, ddc_grid, comm )
  type(t_grid), target, intent(out) :: grid
  character(len=*), intent(in) :: grid_file_name
  type(t_ddc_grid), intent(out), optional :: ddc_grid
  integer, intent(in), optional :: comm
 
  integer :: ierr, fu, d, ie
  integer, allocatable :: t(:,:), e(:,:)
  real(wp), allocatable :: p(:,:)
  character(len=*), parameter :: &
    this_sub_name = 'new_grid_file'
  character(len=1000+len_trim(grid_file_name)) :: message(2)

   ! read the mesh file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(grid_file_name), &
        iostat=ierr,status='old',action='read')
    if(ierr.ne.0) then
      write(message(1),'(a)') 'Problems opening the grid file'
      write(message(2),'(a,a,a)') '  "',trim(grid_file_name),'"'
      call error(this_sub_name,this_mod_name,message(1:2))
    endif
     call read_octave(d,'d',fu)
     call read_octave_al(p,'p',fu) ! allocate( p )
     call read_octave_al(e,'e',fu) ! allocate( e )
     call read_octave_al(t,'t',fu) ! allocate( t )
   close(fu,iostat=ierr)

   call new_grid( grid, d, p, e, t, grid_file_name, &
                  ddc_grid, comm )
   deallocate( p , e , t )

   ! Check if the grid file also contains some additional data, and in
   ! case it does, copy them in the corresponding fields.
   call new_file_unit(fu,ierr)
   open(fu,file=trim(grid_file_name), & ! reopen the grid file
        iostat=ierr,status='old',action='read')
   call locate_var(fu,'edata_legend',ierr)
   if(ierr.eq.0) then ! include the data
     call read_octave_al(grid%edata_legend,'edata_legend',fu)
     call read_octave_al(p,'edata',fu) ! using p as temporary
     d = size(grid%edata_legend)
     do ie=1,grid%ne
       allocate(grid%e(ie)%edata(d))
       grid%e(ie)%edata = p(:,ie)
     enddo
     deallocate( p )
   endif
   close(fu,iostat=ierr)
 
 end subroutine new_grid_file
 
!-----------------------------------------------------------------------

 !> Build a \c t_grid object.
 !!
 !! The main input arguments are
 !! <ul>
 !!  <li> \c d: grid dimension
 !!  <li> \c p: coordinates of the vertices (m,nv)
 !!  <li> \c e: side information (d+3,nb) (see later)
 !!  <li> \c t: element connectivity (d+1,ne) (additional rows ignored)
 !! </ul>
 !! The dimension of the grid is \c d, while the dimension \f$m\f$ of
 !! the space \f$\mathbb{R}^m\f$ where the grid is immersed is
 !! deduced from the number of rows of \c p.
 !!
 !! For compatibility with
 !! <tt>Matlab</tt>\f$^{\textrm{\texttrademark}}\f$, each column of e
 !! has the following structure:
 !! <ul>
 !!  <li> d rows with the vertices of the boundary side
 !!  <li> 2 rows are not used
 !!  <li> boundary marker (scalar integer)
 !!  <li> additional rows are ignored
 !! </ul>
 !!
 !! Concerning the domain decomposition sides, they should not appear
 !! in the \c e array. If there are ddc sides, the related information
 !! is retrieved from the grid file \c grid_file_name (which then must
 !! be present) and returned in the optional argument ddc_grid (which
 !! then must be present as well). If there are no ddc sides, but the
 !! optional argument is present, it is initialized to a consistent,
 !! simplified state, without reading any ddc type information from
 !! the grid file (this makes it possible using this module even for
 !! grids that are not partitioned).
 !<
 subroutine new_grid( grid, d, p, e, t, grid_file_name, &
                      ddc_grid, comm )
  type(t_grid), target, intent(out) :: grid
  integer, intent(in) :: d, t(:,:), e(:,:)
  real(wp), intent(in) :: p(:,:)
  character(len=*), intent(in), optional :: grid_file_name
  type(t_ddc_grid), intent(out), optional :: ddc_grid
  integer, intent(in), optional :: comm

  type t_side_list !< type for the side linked list
    integer, allocatable :: iv(:,:) !< vertices (before reordering)
    type(t_dperm) :: v(2) !< vertices (after reordering)
    integer :: ie(2)      !< elements
    integer :: isl(2)     !< local side indexes
    type(t_side_list), pointer :: next => null()
  end type t_side_list
  type t_e2side !< from element to one side in the side list
    integer, pointer :: ie2   !< opposite element (or boundary marker)
    integer, pointer :: isl2  !< local position on ie2
    integer, pointer :: iv(:) !< vertices (before reordering)
    type(t_perm), pointer :: pi !< permutation to reduced form
  end type t_e2side
  type t_v2se_list !< list for sides/elements connected to a vertex
    integer :: nse = 0
    integer :: ise
    type(t_v2se_list), pointer :: next => null()
  end type t_v2se_list

  integer :: i, ierr, fu, iv, ivl, is, isl, isl2, isb, ie, ie2
  integer :: & ! domain decomposition variables
    ddc_counter, & ! count domain decomposition sides
    ddc_marker,  & ! marker of the ddc boundary (local variable)
    sendb(1), recvb(1) ! buffers must be arrays
  ! heads of the linked list (one list for each vertex)
  type(t_dperm) :: side_list_tail
  type(t_side_list), allocatable, target :: side_list(:)
  type(t_e2side), allocatable :: e2side_list(:,:)
  type(t_v2se_list), allocatable, target :: v2s_list(:), v2e_list(:)
  type(t_ddc_nv), allocatable :: ivl2ivn(:)
  ! error message
  character(len=10000) :: message(5)
  character(len=*), parameter :: &
    this_sub_name = 'new_grid'


   !0) define some simple grid parameters
   ! grid%d has been read directly from the grid file
   grid%d = d
   grid%m = size(p,1)
   call new_me(grid%me,grid%d)
   grid%ne = size(t,2)
   grid%nv = size(p,2)


   ! 1) build the side list (lexicographic order)
   allocate(side_list(grid%nv))
   ! terminate the lists with grid%nv+1, larger than any iv
   side_list_tail = dperm_reduce((/(grid%nv+1, i=1,grid%d)/))
   do iv=1,grid%nv
     allocate(side_list(iv)%next)
     side_list(iv)%next%v(1) = side_list_tail
   enddo
   ! allocate the array used to read the side_list
   allocate(e2side_list(grid%d+1,grid%ne))
   ! build the list: while doing this count the domain decomposition
   ! sides as follows
   ! a) each time a new side is created, consider it as ddc
   ! b) each time a second element is inserted eliminate one ddc
   ! c) for each entry in e, eliminate one ddc
   ! d) count the remaining ddcs
   ddc_counter = 0
   ! At this point we don't know whether there are ddc sides or not,
   ! but we already need the corresponding marker to initialize the
   ! sizes (see insert_side_list). We start considering the following
   ! local value:
   if(size(e,2).gt.0) then ! the boundary could be empty
     ddc_marker = maxval(e(grid%d+3,:)) + 1
   else
     ddc_marker = 1
   endif
   ! The marker is obtained as the first free integer following 
   ! the maximum integer used to denote the domain faces. If
   ! we have 6 faces with 6 different indexes, ddc_marker will be
   ! equal to 7.
   ! Now, we have the following cases:
   ! - if ddc_grid is not present
   !  * if there are no ddc sides, the previous choice for ddc_marker
   !    gives no problems
   !  * if there are some ddc sides, there will be an error
   ! - if ddc_grid is present, we assume that the communicator is also
   !   present and we set ddc_marker to the proper global value; this
   !   will yield the correct results whether or not there are ddc
   !   sides
   if(present(ddc_grid)) then
     if(.not.present(comm)) call error(this_sub_name,this_mod_name, &
                'The optional argument "comm" must be present.')
     ! Make sure all subdomains agree about the value of ddc_marker
     sendb(1) = ddc_marker
     call mpi_allreduce(sendb,recvb,1,mpi_integer,mpi_max,comm,ierr)
     ddc_marker = recvb(1)
   endif
   do ie=1,grid%ne ! element loop
     do isl=1,grid%d+1 ! local side loop
       call insert_side_list( ie , isl , t(grid%me%isv(:,isl),ie) )
     enddo
   enddo
   ! add boundary information
   do isb=1,size(e,2)
     if(e(grid%d+3,isb).le.0) then
       write(message(1),'(a,i7,a,i7)') &
         'Wrong marker ',e(grid%d+3,isb),' on boundary side ',isb
       write(message(2),'(a)') &
         'Boundary markers must be positive.'
       call error(this_sub_name,this_mod_name,message(1:2))
     endif
     call insert_side_list_b(e(grid%d+3,isb),e(1:grid%d,isb))
   enddo
   

   ! 2) allocate the grid (we now know how many ddc)
   !> \bug this should be allocate(t_grid(size(p,1))::grid)
   grid%nb = size(e,2) + ddc_counter
   ! total number of sides
   grid%ns = int( 0.5_wp*real((grid%d+1)*grid%ne+grid%nb,wp) )
   grid%ni = grid%ns - grid%nb

   ! allocate the main grid variables
   allocate( grid%v(grid%nv) )
   do iv=1,grid%nv; allocate( grid%v(iv)%x(grid%m) ); enddo

   allocate( grid%s(grid%ns) )
   do is=1,grid%ns
     allocate( grid%s(is)%iv(grid%d) , grid%s(is)%v(grid%d) , &
       grid%s(is)%xb(grid%m) , grid%s(is)%b(grid%m,grid%d-1) ,&
       grid%s(is)%x0(grid%m) )
   enddo

   allocate( grid%e(grid%ne) )
   do ie=1,grid%ne
     allocate( grid%e(ie)%iv(grid%d+1), grid%e(ie)%v(grid%d+1) ,     &
          grid%e(ie)%is(grid%d+1) ,     grid%e(ie)%s(grid%d+1) ,     &
          grid%e(ie)%pi(grid%d+1) ,     grid%e(ie)%ip(grid%d+1) ,    &
          grid%e(ie)%ie(grid%d+1) ,     grid%e(ie)%e(grid%d+1) ,     &
          grid%e(ie)%iel(grid%d+1) ,    grid%e(ie)%xb(grid%m) ,      &
          grid%e(ie)%b(grid%m,grid%d),  grid%e(ie)%bi(grid%d,grid%m),&
          grid%e(ie)%x0(grid%m)   ,     grid%e(ie)%n(grid%m,grid%d+1))
   enddo


   ! 3) copy back into the grid variable: while doing this, define the
   !  side ordering leaving boundary sides after internal ones
   is = 0
   isb = grid%ni
   allocate( v2s_list(grid%nv) , v2e_list(grid%nv) )
   grid%s%d = -1 ! check that all the sides are initialized
   do ie=1,grid%ne ! element loop

     ! build the element
     grid%e(ie)%d = grid%d
     grid%e(ie)%m = grid%m
     grid%e(ie)%i = ie
     do ivl=1,grid%d+1
       grid%e(ie)%iv(ivl) = t(ivl,ie)
       grid%e(ie)%v(ivl)%p => grid%v( grid%e(ie)%iv(ivl) )
       call v2se_list_add(v2e_list,ie,t(ivl,ie))
     enddo

     ! build the sides
     do isl=1,grid%d+1 ! local side loop

       ie2 = e2side_list(isl,ie)%ie2 ! opposite element
       if(ie2.gt.ie) then ! we have to create a new side

         is = is+1
         isl2 = e2side_list(isl,ie)%isl2 ! local position on ie2

         ! vertex information
         ! edge
         do ivl=1,grid%d ! loop over side vertices
           iv = e2side_list(isl,ie)%iv(ivl)
           call v2se_list_add(v2s_list,is,iv)
         enddo

         ! side information
         grid%s(is)%d = grid%d
         grid%s(is)%m = grid%m
         grid%s(is)%i = is
         ! vertex
         ! the side gets the ordering induced by e1
         grid%s(is)%iv = e2side_list(isl,ie)%iv
         do ivl=1,grid%d
           grid%s(is)%v(ivl)%p => grid%v(grid%s(is)%iv(ivl))
         enddo
         ! element
         grid%s(is)%ie = (/ ie , ie2 /)
         grid%s(is)%e(1)%p => grid%e(ie)
         grid%s(is)%e(2)%p => grid%e(ie2)
         grid%s(is)%isl = (/ isl , isl2 /)

         ! element information
         ! edge
         grid%e(ie)%is(isl) = is
         grid%e(ie)%s(isl)%p => grid%s(is)
         grid%e(ie)%pi(isl) = 1 ! identity
         grid%e(ie)%ip(isl) = &
           idx( grid%me%pi_tab(grid%e(ie)%pi(isl))**(-1) )
         grid%e(ie2)%is(isl2) = is
         grid%e(ie2)%s(isl2)%p => grid%s(is)
         grid%e(ie2)%pi(isl2) = idx(                                 &
           ((e2side_list(isl,ie)%pi)**(-1)) * e2side_list(isl2,ie2)%pi )
         grid%e(ie2)%ip(isl2) = &
           idx( grid%me%pi_tab(grid%e(ie2)%pi(isl2))**(-1) )
         ! neighboring elements
         grid%e(ie)%ie(isl) = ie2
         grid%e(ie)%e(isl)%p => grid%e(ie2)
         grid%e(ie)%iel(isl) = isl2
         grid%e(ie2)%ie(isl2) = ie
         grid%e(ie2)%e(isl2)%p => grid%e(ie)
         grid%e(ie2)%iel(isl2) = isl
 
       elseif(ie2.lt.0) then ! we have to create a boundary side
         
         ! for boundary edges the ordering is provided by e
         isb = isb + 1

         ! vertex information
         ! edge
         do ivl=1,grid%d ! loop over side vertices
           iv = e2side_list(isl,ie)%iv(ivl)
           call v2se_list_add(v2s_list,isb,iv)
         enddo

         ! side information
         grid%s(isb)%d = grid%d
         grid%s(isb)%m = grid%m
         grid%s(isb)%i = isb
         ! vertex
         ! the side gets the ordering induced by e1
         grid%s(isb)%iv = e2side_list(isl,ie)%iv
         do ivl=1,grid%d
           grid%s(isb)%v(ivl)%p => grid%v(grid%s(isb)%iv(ivl))
         enddo
         ! element
         grid%s(isb)%ie = (/ ie , ie2 /) ! ie2 is -marker
         grid%s(isb)%e(1)%p => grid%e(ie)
         grid%s(isb)%e(2)%p => null()
         grid%s(isb)%isl = (/ isl , 0 /)

         ! element information
         ! edge
         grid%e(ie)%is(isl) = isb
         grid%e(ie)%s(isl)%p => grid%s(isb)
         grid%e(ie)%pi(isl) = 1
         grid%e(ie)%ip(isl) = &
           idx( grid%me%pi_tab(grid%e(ie)%pi(isl))**(-1) )
         ! neighboring elements
         grid%e(ie)%ie(isl) = ie2
         grid%e(ie)%e(isl)%p => null()
         grid%e(ie)%iel(isl) = 0

       endif
     enddo
   enddo
   if(any(grid%s%d.lt.0)) then ! something wrong
     write(message(1),'(a)') &
       'The following sides do not belong to any element: '
     i = 1
     do is=1,grid%ns
       if(grid%s(is)%d.lt.0) then
         write(message(2)(i:i+6),'(i7)') is
         i = i+8
         if(i+8.gt.len(message(2))) exit
       endif
     enddo
     write(message(3),'(a)') 'Possible reasons are:'
     write(message(4),'(a)') '  * too many columns in "e"'
     write(message(5),'(a)') '  * zero boundary markers in "e"'
     call error(this_sub_name,this_mod_name,message(1:5))
   endif

   ! Almost done: vertex information
   do iv=1,grid%nv
     grid%v(iv)%m = grid%m
     grid%v(iv)%i = iv
     ! The coordinates should be set in the geometry computations, but
     ! storing them now will allow us to deallocate p together with e
     ! and t, so that it's convenient.
     grid%v(iv)%x = p(:,iv)
     call copy_v2se(iv) ! also set the pointers s and clean the list
   enddo


   ! 4) Domain decomposition optional output
   if(ddc_counter.gt.0) then
     if( (.not.present(ddc_grid)).or. &
         (.not.present(grid_file_name)) ) then
       write(message(1),'(a,a,a,i9,a)') &
         'The grid in "',trim(grid_file_name),'" has ',ddc_counter, &
         ' domain decomposition sides.'
       write(message(2),'(a,a)') &
         ' You must use the optional arguments to retrive the related', &
         ' information.'
       call error(this_sub_name,this_mod_name,message(1:2))
     endif
     call new_file_unit(fu,ierr)
     open(fu,file=trim(grid_file_name), &
          iostat=ierr,status='old',action='read')
       call read_octave_al(ivl2ivn,'ivl2ivn',fu) ! allocate( ivl2ivn )
     close(fu,iostat=ierr)

     ! We can now define the ddc grid
     call compute_ddc( ddc_grid , grid , comm , ddc_marker , ivl2ivn )
   else
     if(present(ddc_grid)) &
       call compute_ddc_simple( ddc_grid , grid , ddc_marker )
   endif


   ! Clean up
   deallocate( e2side_list , v2s_list , v2e_list )
   call side_list_clean()


   ! Compute geometrical quantities
   call compute_geometry(grid)

   ! edata fields are included in  new_grid_file

  contains

   subroutine insert_side_list(ie,isl,isv)
    integer, intent(in) :: ie, isl
    integer, intent(in) :: isv(:)

    type(t_dperm) :: isv_red
    type(t_side_list), pointer :: p, new_side
     
     isv_red = dperm_reduce(isv) ! reduce to basic form
     ! insert in the linked list
     p => side_list(isv_red%x(1)) ! the first vertex indexes the array
     search: do
       ! Notice: the two values v(1) and v(2) have the same field x
       ! (the one used to sort the sides) and different fields pi, so
       ! that the test on the ordering can be done on v(1).
       if(isv_red.lt.p%next%v(1)) then ! add new side to the list
         allocate(new_side)
         allocate(new_side%iv(grid%d,2)); new_side%iv(:,1) = isv
         new_side%v(1)   = isv_red
         new_side%ie(1)  = ie
         new_side%isl(1) = isl
         ! Setting ie(2) as -ddc_marker indicates that the newly
         ! created side is a ddc side. If there will be a second
         ! element, or the side appears in e, this value will be
         ! overwritten, so that after building the side list the ddc
         ! sides can be identified testing ie(2).eq.-ddc_marker
         new_side%ie(2) = -ddc_marker ! default initialization
         e2side_list(isl,ie)%ie2  => new_side%ie(2)
         e2side_list(isl,ie)%isl2 => new_side%isl(2)
         e2side_list(isl,ie)%iv   => new_side%iv(:,1)
         e2side_list(isl,ie)%pi   => new_side%v(1)%pi
         new_side%next   => p%next
         p%next => new_side
         ddc_counter = ddc_counter + 1
         exit search
       elseif(isv_red.eq.p%next%v(1)) then ! add el. to an existing side
         p%next%iv(:,2) = isv
         p%next%v(2)   = isv_red
         p%next%ie(2)  = ie
         p%next%isl(2) = isl
         e2side_list(isl,ie)%ie2  => p%next%ie(1)
         e2side_list(isl,ie)%isl2 => p%next%isl(1)
         e2side_list(isl,ie)%iv   => p%next%iv(:,2)
         e2side_list(isl,ie)%pi   => p%next%v(2)%pi
         ddc_counter = ddc_counter - 1
         exit search
       endif

       p => p%next
     enddo search

   end subroutine insert_side_list

   subroutine insert_side_list_b(marker,isv)
   ! Analogous to insert_side_list, but here we can be sure that the
   ! side is already present.
    integer, intent(in) :: marker, isv(:)

    type(t_dperm) :: isv_red
    type(t_side_list), pointer :: p
     
     isv_red = dperm_reduce(isv) ! reduce to basic form
     ! insert in the linked list
     p => side_list(isv_red%x(1))
     search: do
       if(isv_red.eq.p%next%v(1)) then ! add boundary
         ! check that the same boundary side is not already defined
         error_if: if(p%next%ie(2).ne.-ddc_marker) then
           ! For the check, see the comments in insert_side_list
           write(message(1),'(a)') &
             "Conflicting definition of a boundary side."
           write(message(3),'(a,i5,a)') '(a,',size(isv),'i9)'
           write(message(2),message(3)) &
             "  side vertexes: ", isv
           write(message(3),'(a,i9)') &
             "  ie(2) is already: ", p%next%ie(2)
           call error(this_sub_name,this_mod_name,message(1:3))
         endif error_if
         ! p%next%v(2) left undefined
         p%next%ie(2) = -marker
         ! p%next%isl(2) left undefined
         ddc_counter = ddc_counter - 1
         exit search
       endif

       p => p%next
     enddo search

   end subroutine insert_side_list_b

   subroutine side_list_clean()
    type(t_side_list), pointer :: p, p2

     do iv=1,grid%nv
       p => side_list(iv)%next
       do
        if(.not.associated(p)) exit
         p2 => p
         p => p%next
         deallocate(p2)
       enddo
       side_list(iv)%next => null()
     enddo
     deallocate(side_list)
   end subroutine side_list_clean

   subroutine v2se_list_add(list,ise,iv)
   ! add side/element ise to the list of vertex iv
    type(t_v2se_list), target :: list(:)
    integer, intent(in) :: ise, iv

    type(t_v2se_list), pointer :: p

     p => list(iv)
     do
      if(.not.associated(p%next)) exit
       p => p%next
     enddo
     list(iv)%nse = list(iv)%nse + 1
     allocate(p%next)
     p%next%ise = ise

   end subroutine v2se_list_add

   subroutine copy_v2se(iv)
    integer, intent(in) :: iv

    type(t_v2se_list), pointer :: p, p2

     grid%v(iv)%ns = v2s_list(iv)%nse
     grid%v(iv)%ne = v2e_list(iv)%nse
     allocate( grid%v(iv)%is(grid%v(iv)%ns) , &
               grid%v(iv)%s (grid%v(iv)%ns)  )
     allocate( grid%v(iv)%ie(grid%v(iv)%ne) , &
               grid%v(iv)%e (grid%v(iv)%ne)  )

     p => v2s_list(iv)%next
     do
      if(.not.associated(p)) exit
       grid%v(iv)%is(v2s_list(iv)%nse) = p%ise
       grid%v(iv)%s(v2s_list(iv)%nse)%p => grid%s(p%ise)
       v2s_list(iv)%nse = v2s_list(iv)%nse-1
       p2 => p
       p => p%next
       deallocate(p2)
     enddo

     p => v2e_list(iv)%next
     do
      if(.not.associated(p)) exit
       grid%v(iv)%ie(v2e_list(iv)%nse) = p%ise
       grid%v(iv)%e(v2e_list(iv)%nse)%p => grid%e(p%ise)
       v2e_list(iv)%nse = v2e_list(iv)%nse-1
       p2 => p
       p => p%next
       deallocate(p2)
     enddo

     v2s_list(iv)%next => null()
     v2e_list(iv)%next => null()
   end subroutine copy_v2se

 end subroutine new_grid

!-----------------------------------------------------------------------

 pure subroutine clear_grid( grid )
  type(t_grid), intent(out) :: grid
   call clear( grid%me )
 end subroutine clear_grid
 
!-----------------------------------------------------------------------

 pure subroutine clear_ddc_grid( ddc_grid )
  type(t_ddc_grid), intent(out) :: ddc_grid
 end subroutine clear_ddc_grid
 
!-----------------------------------------------------------------------
 
 !> Affine transformation \f$x = B^E \xi + v^E_0\f$.
 pure function affmap(e,csi) result(x)
  type(t_e), intent(in) :: e
  real(wp), intent(in) :: csi(:,:)
  real(wp) :: x(e%m,size(csi,2))

  integer :: i
 
   x = matmul(e%b,csi)
   do i=1,e%m
     x(i,:) = x(i,:) + e%x0(i)
   enddo

 end function affmap
 
!-----------------------------------------------------------------------

 !> Inverse affine transformation \f$\xi = (B^E)^{-1}(x - v^E_0)\f$.
 pure function iaffmap(e,x) result(csi)
  type(t_e), intent(in) :: e
  real(wp), intent(in) :: x(:,:)
  real(wp) :: csi(e%d,size(x,2))

  integer :: i
  real(wp) :: xx(size(x,1),size(x,2))

   do i=1,e%m
     xx(i,:) = x(i,:) - e%x0(i)
   enddo
   csi = matmul(e%bi,xx)
 
 end function iaffmap
 
!-----------------------------------------------------------------------
 
 !> Compute the diameter \f$h_K\f$ of an element \f$K\f$. This is
 !! defined in terms of the diameters of the sides as \f$h_K =
 !! max_{e\in\partial K}h_e\f$. The case \f$d=1\f$ is handled
 !! separately setting \f$h=|K|\f$.
 !<
 elemental function diameter_e(e) result(h)
  type(t_e), intent(in) :: e
  real(wp) :: h
 
  integer :: i

   if(e%d.eq.1) then
     h = e%vol
   else
     h = 0.0_wp
     do i=1,e%d+1
       h = max( h , diameter(e%s(i)%p) )
     enddo
   endif

 end function diameter_e
 
!-----------------------------------------------------------------------
 
 !> Affine transformation \f$x = B^S \xi + v^S_0\f$.
 pure function affmap_s(s,csi) result(x)
  type(t_s), intent(in) :: s
  real(wp), intent(in) :: csi(:,:)
  real(wp) :: x(s%m,size(csi,2))
 
  integer :: i

   x = matmul(s%b,csi)
   do i=1,s%m
     x(i,:) = x(i,:) + s%x0(i)
   enddo

 end function affmap_s
 
!-----------------------------------------------------------------------

 !> Compute the diameter \f$h_e\f$ of a side. In fact, we use an
 !! approximation, corresponding to assuming that the side is a
 !! regular simplex. This yelds
 !! \f{equation}{
 !!   h_e = \sqrt[d-1]{(d-1)!\,\sqrt{\frac{2^{d-1}}{d}}\,|e|}
 !! \f}
 !! This approximation results in much faster code, and is acceptable
 !! as far as the grid is regular. The case \f$d=1\f$ is handled
 !! separately setting \f$h_e=1\f$.
 !<
 elemental function diameter_s(s) result(h)
  type(t_s), intent(in) :: s
  real(wp) :: h

  integer :: d

   d = s%d
   if(d.eq.1) then
     h = 1.0_wp
   else
     h = (real(fact(d-1),wp)*sqrt(real(2**(d-1),wp)/real(d,wp))*s%a) &
         **(1.0_wp/real(d-1,wp))
   endif

 end function diameter_s
 
!-----------------------------------------------------------------------
 
 !> Linear element size
 !!
 !! Compute a linear dimension characteristic for the element and
 !! proportional to its volume.
 !<
 elemental function el_linear_size(e) result(h)
  type(t_e), intent(in) :: e
  real(wp) :: h

  integer :: d

   d = e%d
   if(d.eq.1) then
     h = 1.0_wp
   else
     h = (real(fact(d),wp)*sqrt(real(2**(d),wp)/real(d+1,wp))*e%vol) &
         **(1.0_wp/real(d,wp))
   endif
 
 end function el_linear_size
 
!-----------------------------------------------------------------------

 !> Element Courant number
 !!
 !! Compute the element Courant number, given a velocity field.
 !!
 !! In general terms, given a velocity vector \f$v\f$, the
 !! corresponding Courant number is defined as
 !! \f{displaymath}{
 !!  C = \frac{\Delta t\,|v|}{h}\qquad{\rm or}\qquad
 !!  C = \frac{\Delta t\,\Phi_v}{|K|},
 !! \f}
 !! where \f$\Delta t\f$ is a time-step, \f$h\f$ is some linear scale
 !! characterizing the element and \f$\Phi_v\f$ is some flux of
 !! \f$v\f$ through the element \f$K\f$. The definition is somewhat
 !! arbitrary, however it should take into account two effects:
 !! <ul>
 !!  <li> the Courant number should depend on the relative orientation
 !!  of \f$v\f$ and \f$K\f$, especially for stretched elements
 !!  <li> sometimes, signal propagation happens in all directions,
 !!  thus one is interested in the maximum value of \f$C\f$ on all
 !!  possible directions of \f$v\f$.
 !! </ul>
 !! One possibility to give a precise definition respecting the above
 !! criteria uses the flux-based definition of the Courant number and
 !! proceeds as follows. 
 !! 
 !! First of all, define the total flux
 !! \f{displaymath}{
 !!  \Phi = \int_{\partial K} v\cdot n\,d\sigma = 
 !!  \int_{\partial K^+} v\cdot n\,d\sigma +
 !!  \int_{\partial K^-} v\cdot n\,d\sigma,
 !! \f}
 !! having introduced the partition of \f$\partial K\f$ into inflow
 !! and outflow regions. Now, observing that the total flux must
 !! vanish, we can define
 !! \f{displaymath}{
 !!  \Phi^+ = 
 !!  \int_{\partial K^+} v\cdot n\,d\sigma =
 !!  - \int_{\partial K^-} v\cdot n\,d\sigma =
 !!  \frac{1}{2}\int_{\partial K} |v\cdot n|\,d\sigma.
 !! \f}
 !! The Courant number is then defined as
 !! \f{displaymath}{
 !!  C = \frac{\Delta t\, \Phi^+}{|K|}.
 !! \f}
 !! \note In this function, the time step is not considered, so the
 !! result is simply \f$\frac{\Phi^+}{|K|}\f$; this quantity should
 !! then be multiplied by \f$\Delta t\f$.
 !!
 !! In fact, using the fact that the sides are straight, we can
 !! compute the flux integrals, obtaining
 !! \f{displaymath}{
 !!  \Phi^+ = \frac{1}{2}\sum_{j=1}^{d+1} |e_j| |n_j\cdot v|.
 !! \f}
 !! Introducing the matrix
 !! \f{displaymath}{
 !!  \Sigma = \left[ \begin{array}{c}
 !!   |e_1|n_1^T \\
 !!   \ldots \\
 !!   |e_{d+1}|n_{d+1}^T \end{array}\right]
 !! \f}
 !! we have
 !! \f{displaymath}{
 !!  \Phi^+ = \frac{1}{2}\| \Sigma\,v\|_1.
 !! \f}
 !! From this expression, we see that the maximum value of
 !! \f$\Phi^+\f$ varying the direction of \f$v\f$ is
 !! \f{displaymath}{
 !!  \Phi^{+,max} = \frac{1}{2}\| \Sigma\|_1 \|v\|_1 = \frac{1}{2}
 !!  \left[ \max_{1\leq j\leq d}\sum_{i=1}^{d+1} |e_i||n_{ij}|\right]
 !!  \left[ \sum_{j=1}^d|v_i|\right].
 !! \f}
 !! Using the fact that
 !! \f{displaymath}{
 !!   |e_1|n_1 = -\sum_{j=2}^{d+1} |e_j|n_j,
 !! \f}
 !! we can also introduce an equivalent quantity 
 !! \f{displaymath}{
 !!  \tilde{\Phi}^+ = \frac{1}{2}\| \tilde{\Sigma}\,v\|_1,
 !! \f}
 !! where
 !! \f{displaymath}{
 !!  \tilde{\Sigma} = \left[ \begin{array}{c}
 !!   |e_2|n_2^T \\
 !!   \ldots \\
 !!   |e_{d+1}|n_{d+1}^T \end{array}\right],
 !! \f}
 !! so that
 !! \f{displaymath}{
 !!  \tilde{\Phi}^+ \leq \Phi^+ \leq 2\tilde{\Phi}^+.
 !! \f}
 !! Now, since
 !! \f{displaymath}{
 !!  |e_j|n_j = det(B)|\hat{e}|B^{-T}\hat{n}_j,
 !! \f}
 !! we see that the matrix \f$\tilde{\Sigma}\f$ is simply
 !! \f$-det(B)|\hat{e}|B^{-1}\f$, and using the equivalence between
 !! \f$\|\cdot\|_1\f$ and
 !! \f$\|\cdot\|_2\f$ we can define a third equivalent flux
 !! \f{displaymath}{
 !!  \tilde{\tilde{\Phi}}^+ = \frac{|det(B)||\hat{e}|}{2}\|
 !!  B^{-1}\,v\|_2.
 !! \f}
 !! This leads to the equivalent Courant number
 !! \f{displaymath}{
 !!  \tilde{\tilde{C}} = \frac{d}{2} \sqrt{ \sum_{j=1}^d \left(
 !!  \frac{\Delta t\,v^K_j}{\lambda^K_j} \right)^2},
 !! \f}
 !! where \f$\lambda^K_j\f$ are the principal deformations of \f$B\f$
 !! and \f$v^K_j\f$ are the components of \f$v\f$ along the
 !! corresponding principal directions, and we have used
 !! \f{displaymath}{
 !!  \frac{|det(B)||\hat{e}|}{2|K|} =
 !!  \frac{|\hat{e}|}{2|\hat{K}|} =
 !!  \frac{d!}{2 (d-1)!} = \frac{d}{2}.
 !! \f}
 !! Notice that this definition is in fact based on a set of
 !! characteristic element lengths, namely the \f$\lambda^K_j\f$.
 !!
 !! Summarizing, this function computes three Courant numbers:
 !! \f{displaymath}{
 !!  C = \frac{1}{2|K|}\|\Sigma\,v\|_1, \qquad
 !!  \tilde{\tilde{C}} = \frac{d}{2}\|B^{-1}\,v\|_2, \qquad
 !!  C^{max} = \frac{1}{2|K|}\|\Sigma\|_1\|\,v\|_1.
 !! \f}
 !! The fourth Courant number \f$\tilde{\tilde{C}}^{max}\f$ is not
 !! computed, since it would require \f$\|B^{-1}\|_2\f$.
 pure subroutine el_courant(c,ct,cm,e,esa,v,a)
  type(t_e), intent(in) :: e !< element
  !> side areas \f$|e_j|\f$; this argument is introduced here to avoid
  !! passing the whole grid, it can be built on-the-fly as
  !! \code
  !!  esa = (/ ( e%s(i)%p%a  , i=1,e%d+1 ) /)
  !! \endcode
  real(wp), intent(in) :: esa(:)
  real(wp), intent(in), optional :: &
    v(:), & !< used to compute \c c and \c ct
    a       !< used to compute \f$C^{max}\f$, with \f$\|v\|_1=a\f$
  real(wp), intent(out), optional :: &
    c,  & !< \f$C\f$
    ct, & !< \f$\tilde{\tilde{C}}\f$
    cm    !< \f$C^{max}\f$

   if(present(v)) then
     c  = (0.5_wp/e%vol) * sum(abs(esa*matmul(v,e%n)))
     ct = (0.5_wp*real(e%d,wp)) * sqrt(sum(matmul(e%bi,v)**2))
   endif

   if(present(a)) &
     cm = (0.5_wp/e%vol) * maxval(matmul(abs(e%n),esa))*abs(a)

 end subroutine el_courant
 
!-----------------------------------------------------------------------
 
 !> Use a search algorithm to locate the element a point belongs to
 !!
 !! The indexes of the points that are found within the given grid are
 !! returned in \c ip; hence there are no problems if a requested
 !! point does not belong to the given grid (useful when working with
 !! decomposed grids).
 !<
 pure subroutine locate_point(ie,x,grid,ie0,ip)
  integer, allocatable, intent(out) :: ie(:) !< elements
  real(wp),     intent(in) :: x(:,:) !< points as column vectors
  type(t_grid), intent(in) :: grid   !< grid
  integer,              intent(in), optional :: ie0(:) !< initial guess
  integer, allocatable, intent(out), optional :: ip(:) !< loc. points
 
  real(wp), parameter :: toll = 100.0_wp*epsilon(0.0_wp)
  logical :: search_completed, need_restart, found_restart
  logical, allocatable :: checked(:)
  integer :: np, i, je
  integer, allocatable :: wip(:), wie(:)
  real(wp) :: xi(grid%d+1)
  character(len=*), parameter :: &
    this_sub_name = 'locate_point'

   np = size(x,2)
   allocate( wip(np) , wie(np) ) ! work arrays
   allocate( checked(grid%ne) )

   np = 0 ! counter
   point_do: do i=1,size(x,2)

     checked = .false.

     if(present(ie0)) then
       je = ie0(i)
     else
       je = 1
     endif
     search_completed = .false.
     do
      if(search_completed) exit

       ! barycentric coordinates (see also iaffmap)
       xi(2:) = matmul( grid%e(je)%bi , x(:,i)-grid%e(je)%x0 )
       xi(1)  = 1.0_wp - sum(xi(2:))
       if(minval(xi).ge.-toll) then ! found the element
         np = np+1
         wip(np) = i
         wie(np) = je
         search_completed = .true.
       else
         checked(je) = .true.
         je = grid%e(je)%ie( minloc(xi,dim=1) )
         if(je.lt.1) then
           need_restart = .true.
         elseif(checked(je)) then
           need_restart = .true.
         else
           need_restart = .false.
         endif
         if(need_restart) then
           call get_next(je,found_restart)
           if(.not.found_restart) then ! the point can not be located
             search_completed = .true.
           endif
         endif
       endif
     enddo
     
   enddo point_do

   allocate( ie(np) )
   ie = wie(1:np)
   if(present(ip)) then
     allocate( ip(np) )
     ip = wip(1:np)
   endif

   deallocate( wip , wie , checked )
 contains

  pure subroutine get_next(ie,found)
   integer, intent(out) :: ie
   logical, intent(out) :: found

    found = .false.
    do ie=1,grid%ne
      if(.not.checked(ie)) then
        found = .true.
        exit
      endif
    enddo
  end subroutine get_next

 end subroutine locate_point
 
!-----------------------------------------------------------------------
 
 !> Set the domain decomposition data structures
 subroutine compute_ddc( ddc_grid , grid , comm , ddc_marker , ivl2ivn )
  type(t_ddc_grid), intent(out) :: ddc_grid
  type(t_grid), intent(in) :: grid
  integer, intent(in) :: comm, ddc_marker
  type(t_ddc_nv), intent(in) :: ivl2ivn(:)

  integer :: i, id, iv, ivl, ivn, ivseg(3), is, isn, ie, dim1, ierr
  integer :: req
  integer, allocatable :: sendbuf(:,:), recvbuf(:,:), dim2_id(:), &
    is_s(:), reqs(:), mpi_stat(:,:), ids(:)
  type(t_dperm) :: p1, p2
  character(len=1000) :: message(1)
  character(len=*), parameter :: &
    this_sub_name = 'compute_ddc'

   ! 1) General setup: nd, id, ddc_marker
   call mpi_comm_size(comm,ddc_grid%nd,ierr)
   call mpi_comm_rank(comm,ddc_grid%id,ierr)
   allocate( ddc_grid%nnv_id(0:ddc_grid%nd-1) , &
             ddc_grid%nns_id(0:ddc_grid%nd-1) )
   ddc_grid%nnv_id = 0
   ddc_grid%nns_id = 0
   ddc_grid%ddc_marker = ddc_marker

   ! 2) Global entity information
   allocate( ddc_grid%gv(grid%nv) , &
             ddc_grid%gs(grid%ns) , &
             ddc_grid%ge(grid%ne) )

   ! 3) Neighbour entity information: nnv, nns, gv%ni, gs%ni
   ddc_grid%nnv = 0
   ddc_grid%nns = 0
   do is=grid%ni+1,grid%ns
     if(grid%s(is)%ie(2).eq.-ddc_grid%ddc_marker) then
       ddc_grid%nns = ddc_grid%nns + 1
       ddc_grid%gs(is)%ni = ddc_grid%nns
       do ivl=1,grid%d
         iv = grid%s(is)%iv(ivl)
         if(ddc_grid%gv(iv)%ni.eq.0) then
           ddc_grid%nnv = ddc_grid%nnv + 1
           ddc_grid%gv(iv)%ni = ddc_grid%nnv
         endif
       enddo
     endif
   enddo
   allocate( ddc_grid%nv(ddc_grid%nnv) , &
             ddc_grid%ns(ddc_grid%nns) )

   ! 4) Neighbour vertexes: nv%all_fields, nnv_id
   do ivn=1,size(ivl2ivn)
     ! This works in both cases: whether ivl2ivn contains only ddc
     ! nodes or all the nodes (which depends on the input file)
     if(ivl2ivn(ivn)%nd.gt.0) then
       iv = ddc_grid%gv(ivl2ivn(ivn)%i)%ni
       ddc_grid%nv(iv) = ivl2ivn(ivn)
       ! recall that nv%id is an array
       ddc_grid%nnv_id(ddc_grid%nv(iv)%id) =   &
         ddc_grid%nnv_id(ddc_grid%nv(iv)%id) + 1
     endif
   enddo
   ! Check that the input file was complete
   do iv=1,grid%nv
     if(ddc_grid%gv(iv)%ni.gt.0) then
       ivn = ddc_grid%gv(iv)%ni
       if(.not.allocated(ddc_grid%nv(ivn)%id)) then
         write(message(1),'(a,i9,a,i9,a)') &
       'Domain decomposition information is missing for vertex ' , &
       iv,' on subdomain ',ddc_grid%id,'.'
         call error(this_sub_name,this_mod_name,message(1))
       endif
     endif
   enddo

   ! 5) Neighbour Side information

   ! 5.1) Neighbour sides: ns%i, ns%is, nns_id
   call ns_disambiguate('init',grid,ddc_grid)
   do is=grid%ni+1,grid%ns
     if(ddc_grid%gs(is)%ni.ne.0) then
       ddc_grid%ns(ddc_grid%gs(is)%ni)%i  = is
       ! The id field is initialized to -1 to check later whether a
       ! side has more that one neighbour, which indicates an error in
       ! the grid files.
       ddc_grid%ns(ddc_grid%gs(is)%ni)%id = -1

       call get_siden_id(ids,ddc_grid%nv(ddc_grid%gv(grid%s(is)%iv)%ni))
       if(size(ids).eq.1) then
         id = ids(1) ! usually this is the case
         ddc_grid%ns(ddc_grid%gs(is)%ni)%id = id
         ! the remaining fields require a communication
         ddc_grid%nns_id(id) = ddc_grid%nns_id(id)+1
       else
         if(size(ids).eq.0) then
           write(message(1),*) "Can't find a neighbour for side ",is
           call error(this_sub_name,this_mod_name,message(1))
         else ! we need to disambiguate
           call ns_disambiguate('insert',grid,ddc_grid,is=is,ids=ids)
         endif
       endif
     endif
   enddo
   call ns_disambiguate('disambiguate',grid,ddc_grid,comm=comm)

   ! 5.2) all to all communication: ns%in, ns%p_s2s
   dim1 = grid%d+1
   allocate(sendbuf(dim1,ddc_grid%nns))
   allocate(recvbuf(dim1,ddc_grid%nns))
   allocate(is_s(0:ddc_grid%nd-1))
   do i=0,ddc_grid%nd-1 ! set the offset array
     is_s(i) = sum(ddc_grid%nns_id(0:i-1))
   enddo
   do isn=1,ddc_grid%nns
     id = ddc_grid%ns(isn)%id
     is_s(id) = is_s(id) + 1
     ! first lines: neighbour vertex indexes
     do ivl=1,grid%d
       ivn = ddc_grid%gv(grid%s(ddc_grid%ns(isn)%i)%iv(ivl))%ni
       i = findloc( ddc_grid%nv(ivn)%id , id , dim=1 )
       sendbuf(ivl,is_s(id)) = ddc_grid%nv(ivn)%in(i)
     enddo
     ! last line: local side index
     sendbuf(dim1 ,is_s(id)) = ddc_grid%ns(isn)%i
   enddo
   do i=0,ddc_grid%nd-1 ! reset the offset array
     is_s(i) = sum(ddc_grid%nns_id(0:i-1))
   enddo
   call mpi_alltoallv(                                           &
          sendbuf, dim1*ddc_grid%nns_id, dim1*is_s, mpi_integer, &
          recvbuf, dim1*ddc_grid%nns_id, dim1*is_s, mpi_integer, &
                      comm,ierr)
   do isn=1,ddc_grid%nns ! loop on the received buffer
     ! local side index
     is = get_siden_is(grid%v(recvbuf(1:grid%d,isn)))
     ddc_grid%ns(ddc_grid%gs(is)%ni)%in = recvbuf(grid%d+1,isn)
     ! now find the permutation: s%iv -> v
     p1 = dperm_reduce(grid%s(is)%iv)
     p2 = dperm_reduce(recvbuf(1:grid%d,isn))
     ddc_grid%ns(ddc_grid%gs(is)%ni)%p_s2s = idx( (p2%pi**(-1))*p1%pi )
   enddo
   deallocate(sendbuf,recvbuf,is_s)

   ! 6) Global vertex/side/element numbering: gv%gi, gs%gi, ge%gi

   ! This is done in three steps:
   ! 1) get the indexes from the subdomains with lower id 
   ! 2) assign the remaining indexes
   ! 3) send the indexes to the subdomains with higher id

   ! To allow a correct global numbering of the degrees of freedom,
   ! each domain must also tell the next one what is its maximum
   ! global index. This is done with the following send/recv.
   if(ddc_grid%id.gt.0) then
     call mpi_irecv( ivseg(1) , 3 , mpi_integer ,         &
                     ddc_grid%id-1, 1 , comm , req , ierr )
   else
     ivseg = 0
   endif
   allocate(dim2_id(0:ddc_grid%id-1))
   dim2_id = ddc_grid%nnv_id(0:ddc_grid%id-1) + &
             ddc_grid%nns_id(0:ddc_grid%id-1) 
   allocate( reqs    (                0:ddc_grid%id-1) , &
             mpi_stat(mpi_status_size,0:ddc_grid%id-1) )
   allocate(recvbuf(2,sum(dim2_id)))
   reqs = mpi_request_null
   is = 1
   do i=0,ddc_grid%id-1
     if(dim2_id(i).gt.0) &
       call mpi_irecv( recvbuf(1,is) , 2*dim2_id(i) , mpi_integer , &
                       i, 2, comm, reqs(i), ierr )
     is = is + dim2_id(i)
   enddo
   if(ddc_grid%id.gt.0) call mpi_wait( req , mpi_stat(:,0) , ierr )
   call mpi_waitall( size(reqs) , reqs , mpi_stat , ierr )
   deallocate(reqs,mpi_stat)
   ddc_grid%gv%gi = 0
   ddc_grid%gs%gi = 0
   ddc_grid%ge%gi = 0
   allocate(is_s(0:ddc_grid%id-1))
   do i=0,ddc_grid%id-1 ! set the vertex offsets
     is_s(i) = sum( dim2_id(:i-1) )
   enddo
   do i=0,ddc_grid%id-1
     do is=is_s(i)+1,is_s(i)+ddc_grid%nnv_id(i)
       iv = recvbuf(1,is)
       if(iv.gt.grid%nv) then ! somebody is sending wrong data
         write(message(1),'(a,i10,a,i10,a,i10,a,a,i10)')               &
           'Subdomain ',i,' is sending local index ',iv,               &
           ' (position in the two-procs buffer ',is-is_s(i),')',       &
           ' which is larger than the number of local vertexes ',grid%nv
         call error(this_sub_name,this_mod_name,message)
       endif
       if(ddc_grid%gv(iv)%gi.eq.0) then ! define the global index
         ddc_grid%gv(iv)%gi = recvbuf(2,is)
       else ! in principle, nothing to do, however better to check
         if(ddc_grid%gv(iv)%gi.ne.recvbuf(2,is)) &
           call error(this_sub_name,this_mod_name, &
  'Two neighbours are trying to define two different global indexes.')
       endif
     enddo
   enddo
   ! set the side offsets
   is_s = is_s + ddc_grid%nnv_id(0:ddc_grid%id-1)
   do i=0,ddc_grid%id-1
     do is=is_s(i)+1,is_s(i)+ddc_grid%nns_id(i)
       iv = recvbuf(1,is) ! is already used as shift
       ! This is simpler than for the vertexes, since each side is
       ! shared by exactly two subdomains.
       ddc_grid%gs(iv)%gi = recvbuf(2,is)
     enddo
   enddo
   deallocate(dim2_id,is_s,recvbuf)
   allocate(dim2_id(ddc_grid%id+1:ddc_grid%nd-1))
   dim2_id = ddc_grid%nnv_id(ddc_grid%id+1:ddc_grid%nd-1) + &
             ddc_grid%nns_id(ddc_grid%id+1:ddc_grid%nd-1) 
   allocate(sendbuf(2,sum(dim2_id)))
   allocate(is_s(ddc_grid%id+1:ddc_grid%nd-1))
   do i=ddc_grid%id+1,ddc_grid%nd-1 ! set the offset array
     is_s(i) = sum( dim2_id(:i-1) )
   enddo
   do iv=1,grid%nv
     if(ddc_grid%gv(iv)%gi.eq.0) then ! define the global index
       ivseg(1) = ivseg(1)+1
       ddc_grid%gv(iv)%gi = ivseg(1)
     endif
     if(ddc_grid%gv(iv)%ni.ne.0) then ! possibly fill in sendbuf
       ivn = ddc_grid%gv(iv)%ni
       do i=1,ddc_grid%nv(ivn)%nd
         id = ddc_grid%nv(ivn)%id(i)
         if(id.gt.ddc_grid%id) then ! send this info
           is_s(id) = is_s(id)+1
           sendbuf(1,is_s(id)) = ddc_grid%nv(ddc_grid%gv(iv)%ni)%in(i)
           sendbuf(2,is_s(id)) = ddc_grid%gv(iv)%gi
         endif
       enddo
     endif
   enddo
   do is=1,grid%ns
     if(ddc_grid%gs(is)%gi.eq.0) then ! define the global index
       ivseg(2) = ivseg(2)+1
       ddc_grid%gs(is)%gi = ivseg(2)
     endif
     if(ddc_grid%gs(is)%ni.ne.0) then ! possibly fill in sendbuf
       isn = ddc_grid%gs(is)%ni
       id  = ddc_grid%ns(isn)%id
       if(id.gt.ddc_grid%id) then ! send this info
         is_s(id) = is_s(id)+1
         sendbuf(1,is_s(id)) = ddc_grid%ns(ddc_grid%gs(is)%ni)%in
         sendbuf(2,is_s(id)) = ddc_grid%gs(is)%gi
       endif
     endif
   enddo
   do ie=1,grid%ne
     ivseg(3) = ivseg(3)+1
     ddc_grid%ge(ie)%gi = ivseg(3)
   enddo
   deallocate(is_s)
   allocate( reqs    (                ddc_grid%id+1:ddc_grid%nd-1) , &
             mpi_stat(mpi_status_size,ddc_grid%id+1:ddc_grid%nd-1) )
   if(ddc_grid%id.lt.ddc_grid%nd-1) &
     call mpi_isend( ivseg(1) , 3 , mpi_integer , ddc_grid%id+1 , 1 , &
                     comm , req , ierr )
   reqs = mpi_request_null
   is = 1
   do i=ddc_grid%id+1,ddc_grid%nd-1
     if(dim2_id(i).gt.0) &
       call mpi_isend( sendbuf(1,is) , 2*dim2_id(i) , mpi_integer , &
                      i, 2, comm, reqs(i) , ierr )
     is = is + dim2_id(i)
   enddo
   ! We have to wait before deallocating the send buffer
   deallocate(dim2_id)

   ! 7) Almost done, the last subdomain must update the others on the
   ! total number of generated entities
   allocate(dim2_id(3)) ! used as temporary
   dim2_id(1) = maxval(ddc_grid%gv%gi)
   dim2_id(2) = maxval(ddc_grid%gs%gi)
   dim2_id(3) = maxval(ddc_grid%ge%gi)
   call mpi_bcast(dim2_id,3,mpi_integer,ddc_grid%nd-1,comm,ierr)
   ddc_grid%ngv = dim2_id(1)
   ddc_grid%ngs = dim2_id(2)
   ddc_grid%nge = dim2_id(3)
   deallocate(dim2_id)

   if(ddc_grid%id.lt.ddc_grid%nd-1) &
     call mpi_wait( req , mpi_stat(:,ddc_grid%id+1) , ierr )
   call mpi_waitall( size(reqs) , reqs , mpi_stat , ierr )
   deallocate(sendbuf,reqs,mpi_stat)

 contains

  !> The id of the side neighboring subdomain is obtained as the
  !! unique intersection of the id-s of the neighboring subdomains of
  !! the vertexes.
  !!
  !! It is possible (albeit rare) that more than two domains share all
  !! the vertexes of a face, even if only two of them can share the
  !! face. This subroutine returns the id of all the subdomains wich
  !! share all the face vertexes. Notice that the only way to find out
  !! which subdomain really shares the face is by communication.
  !<
  pure subroutine get_siden_id(id,vid)
   integer, allocatable, intent(out) :: id(:)
   type(t_ddc_nv),       intent(in)  :: vid(:) !< side vertexes

   logical :: found
   integer :: i, ii, n, wid, wids(vid(1)%nd)

    n = 0
    search_do: do i=1,vid(1)%nd
      wid = vid(1)%id(i) ! candidate domain
      found = .true.
      check_do: do ii=2,size(vid)
        found = any(vid(ii)%id.eq.wid)
        if(.not.found) exit check_do
      enddo check_do
      if(found) then
        n = n+1
        wids(n) = wid
      endif
    enddo search_do

    allocate(id(n)); id = wids(1:n)
  end subroutine get_siden_id

  !> The is of the side is obtained as the unique intersection of the
  !! is-s of the sides connected with the specified vertexes.
  !<
  pure function get_siden_is(v) result(is)
   type(t_v), intent(in) :: v(:) !< vertexes
   integer :: is !< side index

   logical :: found
   integer :: i, ii

    search_do: do i=1,v(1)%ns
      is = v(1)%is(i)
      found = .true.
      check_do: do ii=2,size(v)
        found = found .and. any(v(ii)%is.eq.is)
        if(.not.found) exit check_do
      enddo check_do
      if(found) exit search_do
    enddo search_do

    if(.not.found) is = huge(is) ! error
  end function get_siden_is

  pure function findloc(v,x,dim) result(i)
  ! This function should be eliminated once the new intrinsic with the
  ! same name is supported
   integer, intent(in) :: v(:), x
   integer, intent(in), optional :: dim
   integer :: i
   integer :: j
    i = 0
    search_do: do j=1,size(v)
      if(v(j).eq.x) then
        i = j
        exit search_do
      endif
    enddo search_do
  end function findloc

 end subroutine compute_ddc

!-----------------------------------------------------------------------

 subroutine compute_ddc_simple( ddc_grid , grid , ddc_marker )
  type(t_ddc_grid), intent(out) :: ddc_grid
  type(t_grid), intent(in) :: grid
  integer, intent(in) :: ddc_marker

  integer :: i

   ddc_grid%nd = 1
   ddc_grid%id = 0
   ddc_grid%ngv = grid%nv
   ddc_grid%ngs = grid%ns
   ddc_grid%nge = grid%ne
   ddc_grid%nnv = 0
   ddc_grid%nns = 0
   allocate(ddc_grid%nnv_id(0:0),ddc_grid%nns_id(0:0))
   ddc_grid%nnv_id = 0; ddc_grid%nns_id = 0
   ddc_grid%ddc_marker = ddc_marker
   allocate(ddc_grid%gv(grid%nv),ddc_grid%gs(grid%ns), &
            ddc_grid%ge(grid%ne))
   ddc_grid%gv%gi = (/(i, i=1,grid%nv)/)
   ddc_grid%gs%gi = (/(i, i=1,grid%ns)/)
   ddc_grid%ge%gi = (/(i, i=1,grid%ne)/)
   allocate(ddc_grid%nv(0),ddc_grid%ns(0))

 end subroutine compute_ddc_simple

!-----------------------------------------------------------------------

 subroutine ns_disambiguate(action,grid,ddc_grid,is,ids,comm)
  character(len=*), intent(in) :: action
  type(t_grid),     intent(in) :: grid
  type(t_ddc_grid), intent(inout) :: ddc_grid
  integer, intent(in), optional :: is, ids(:), comm

  type t_dis
    integer :: counter
    integer, allocatable :: is(:) ! sides to check with a neighbour
  end type t_dis
  type(t_dis), allocatable, save :: dis_status(:)
  character(len=*), parameter :: this_sub_name = "ns_disambiguate"

  integer :: id, i, j, k, inv, ierr
  integer, allocatable :: tmp(:), rec_counter(:), senoff(:), &
    recoff(:), sendbuf(:,:), recvbuf(:,:)
  character(len=1000) :: message(2)

   action_case: select case(action)

    case('init')
     allocate(dis_status(0:ddc_grid%nd-1))
     do id=0,ddc_grid%nd-1
       dis_status(id)%counter = 0
       allocate(dis_status(id)%is(0))
     enddo

    case('insert')
     do i=1,size(ids) ! loop on the candidate neighboring subdomains
       id = ids(i)
       allocate( tmp(dis_status(id)%counter) )
       tmp = dis_status(id)%is
       deallocate(dis_status(id)%is)
       dis_status(id)%counter = dis_status(id)%counter + 1
       allocate(dis_status(id)%is(dis_status(id)%counter))
       dis_status(id)%is(1:dis_status(id)%counter-1) = tmp
       dis_status(id)%is(  dis_status(id)%counter  ) = is
       deallocate(tmp)
     enddo

    case('disambiguate') ! and clean
     ! we don't know how much data we are going to receive
     allocate(rec_counter(0:ddc_grid%nd-1))
     call mpi_alltoall( dis_status%counter , 1 , mpi_integer , &
                        rec_counter        , 1 , mpi_integer , &
                        comm, ierr )
     ! now we can send/receive the side data
     allocate(senoff(0:ddc_grid%nd-1),recoff(0:ddc_grid%nd-1))
     senoff(0) = 0; recoff(0) = 0
     do id=1,ddc_grid%nd-1
       senoff(id) = senoff(id-1) + dis_status(id-1)%counter
       recoff(id) = recoff(id-1) + rec_counter(id-1)
     enddo
     allocate( sendbuf(grid%d,sum(dis_status%counter)) , &
               recvbuf(grid%d,sum(   rec_counter    )) )
     j = 0
     do id=0,ddc_grid%nd-1
       do i=1,dis_status(id)%counter
         j = j+1
         do k=1,grid%d
           inv = ddc_grid%gv(grid%s(dis_status(id)%is(i))%iv(k))%ni
           sendbuf(k,j) = ddc_grid%nv(inv)%in(   &
             findloc( ddc_grid%nv(inv)%id , id ) )
         enddo
       enddo
     enddo
     call mpi_alltoallv(                                               &
       sendbuf, grid%d*dis_status%counter, grid%d*senoff, mpi_integer, &
       recvbuf, grid%d*   rec_counter    , grid%d*recoff, mpi_integer, &
                        comm,ierr)
     ! send/receive the feedback
     deallocate(sendbuf); allocate( sendbuf(1,sum(rec_counter)) )
     do i=1,size(recvbuf,2) ! loop on the received data
       sendbuf(1,i) = check_side_exists(grid%v(recvbuf(:,i)))
     enddo
     deallocate(recvbuf); allocate( recvbuf(1,sum(dis_status%counter)) )
     call mpi_alltoallv(                                               &
       sendbuf,    rec_counter    , recoff, mpi_integer, &
       recvbuf, dis_status%counter, senoff, mpi_integer, &
                        comm,ierr)
     j = 0
     do id=0,ddc_grid%nd-1
       do i=1,dis_status(id)%counter
         j = j+1
         if(recvbuf(1,j).eq.1) then ! positive feedback
           ! Make sure that each side has at most one positive feedback
           k = ddc_grid%ns(ddc_grid%gs(dis_status(id)%is(i))%ni)%id
           if(k.ge.0) then
             write(message(1),'(a,i6)')      &
               'Unable to identify a neighbour on side ', &
               dis_status(id)%is(i)
             write(message(2),'(a,i4,a,i4)') &
               '  possible candidates are (at least) ',k,' and ',id
             call error(this_sub_name,this_mod_name,message(1:2))
           endif
           ddc_grid%ns(ddc_grid%gs(dis_status(id)%is(i))%ni)%id = id
           ! the remaining fields require a communication
           ddc_grid%nns_id(id) = ddc_grid%nns_id(id)+1
         endif
       enddo
     enddo

     deallocate(dis_status)
     deallocate(rec_counter,senoff,recoff,sendbuf,recvbuf)

   end select action_case

 contains

  !> A side is checked as the (unique, if any) intersection of the
  !! sides connected to the d vertexes. Notice that this makes sense
  !! only if called with d vertexes.
  !<
  pure function check_side_exists(v) result(yesno)
   type(t_v), intent(in) :: v(:)
   integer :: yesno

   logical :: found
   integer :: i, is, ii
    
    yesno = 0
    search_do: do i=1,v(1)%ns
      is = v(1)%is(i) ! candidate side
      found = .false.
      check_do: do ii=2,size(v)
        found = any(v(ii)%is.eq.is)
        if(.not.found) exit check_do
      enddo check_do
      if(found) then
        yesno = 1
        exit search_do
      endif
    enddo search_do
  end function check_side_exists

  pure function findloc(v,x,dim) result(i)
  ! This function should be eliminated once the new intrinsic with the
  ! same name is supported
   integer, intent(in) :: v(:), x
   integer, intent(in), optional :: dim
   integer :: i
   integer :: j
    i = 0
    search_do: do j=1,size(v)
      if(v(j).eq.x) then
        i = j
        exit search_do
      endif
    enddo search_do
  end function findloc

 end subroutine ns_disambiguate

!-----------------------------------------------------------------------
 
 !> Compute all the geometric properties of the grid.
 pure subroutine compute_geometry( grid )
  type(t_grid), target, intent(inout) :: grid
 
  integer :: is, ie, i
  integer, allocatable :: pp(:)
  real(wp), allocatable :: qq(:,:), rr(:,:), rri(:,:)
  type(t_s), pointer :: s
  type(t_e), pointer :: e
  character(len=*), parameter :: &
    this_sub_name = 'compute_geometry'

   ! vertices are completed already

   ! element loop
   allocate(     pp(grid%d)     , qq(grid%m,grid%d) , &
             rr(grid%d,grid%d) , rri(grid%d,grid%d) )
   do ie=1,grid%ne
     e => grid%e(ie)

     e%xb = 0.0_wp
     do i=1,e%d+1
       e%xb = e%xb + e%v(i)%p%x
     enddo
     e%xb = e%xb/real(e%d+1,wp)

     ! the matrix b is the main information source for the geometry
     do i=1,e%d
       e%b(:,i) = e%v(i+1)%p%x - e%v(1)%p%x
     enddo
     ! The metric is   G = B^T * B   or, introducing the QR
     ! factorization of B,   B = Q * R,   G = R^T * R. This implies:
     !   det(G) = det(R)^2
     !   B^(-1) = R^(-1) * Q^T
     ! and taking into account the permutation   B * P = Q * R   the
     ! second relation becomes
     !   P * B^(-1) = R^(-1) * Q^T
     qq      = e%b; call qr(qq,rr,pp)
     e%vol   = 1.0_wp/real(fact(e%d),wp) ! det(R) = prod(diag(R))
     do i=1,e%d; e%vol = e%vol*rr(i,i); enddo
     e%det_b = e%vol/grid%me%vol
     call invmat(rr,rri)
     e%bi(pp,:) = matmul(rri,transpose(qq))
     e%x0    = e%v(1)%p%x
     ! Notice: it turns out that the covariant basis (in Cartesian
     ! coordinates) is provided by the rows of e%bi
     e%n(:,1) = sum(e%bi,1) / sqrt(real(e%d,wp))
     do i=2,e%d+1
       e%n(:,i) = - e%bi(i-1,:)
     enddo
     ! normalization
     do i=1,e%d+1
       e%n(:,i) = e%n(:,i) / sqrt(sum(e%n(:,i)**2))
     enddo
   enddo
   deallocate( pp , qq , rr , rri )

   ! side loop
   allocate(     pp(grid%d-1)     ,     qq(grid%m,grid%d-1) , &
             rr(grid%d-1,grid%d-1) )
   do is=1,grid%ns
     s => grid%s(is)

     s%xb = 0.0_wp
     do i=1,s%d
       s%xb = s%xb + s%v(i)%p%x
     enddo
     s%xb = s%xb/real(s%d,wp)

     ! the matrix b is the main information source for the geometry
     do i=1,s%d-1
       s%b(:,i) = s%v(i+1)%p%x - s%v(1)%p%x
     enddo
     ! the metric is   G = B^T * B   and everything proceeds as before
     qq = s%b; call qr(qq,rr,pp)
     s%a = 1.0_wp/real(fact(s%d-1),wp)
     do i=1,s%d-1; s%a = s%a*rr(i,i); enddo
     s%x0    = s%v(1)%p%x
   enddo
   deallocate( pp , qq , rr )

   grid%vol = sum(grid%e%vol)
 
 end subroutine compute_geometry
 
!-----------------------------------------------------------------------

 !> Build a \f$(d-1)\f$-dimensional grid from a collection of faces of
 !! a \f$d\f$-dimensional grid.
 !!
 !! \note The vertex ordering of the elements of the subgrid is
 !! automatically adapted to the order defined by the sides of the
 !! parent grid.
 !<
 subroutine new_subgrid( sub_grid, sub_v2v, sub_e2s, &
                         grid, grid_file_name )
  type(t_grid), target, intent(out) :: sub_grid
  integer, allocatable, intent(out) :: sub_v2v(:), sub_e2s(:)
  type(t_grid), intent(in) :: grid
  character(len=*), intent(in) :: grid_file_name

  logical :: found
  integer :: fu, ierr, nv, ie, iv, isl, is, ib
  integer, allocatable :: e(:,:), t(:,:), v2v(:)
  real(wp), allocatable :: p(:,:)
  character(len=1000+len_trim(grid_file_name)) :: message(2)
  character(len=*), parameter :: &
    this_sub_name = 'new_subgrid'

   ! read the mesh file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(grid_file_name), &
        iostat=ierr,status='old',action='read')
    if(ierr.ne.0) then
      write(message(1),'(a)') 'Problems opening the grid file'
      write(message(2),'(a,a,a)') '  "',trim(grid_file_name),'"'
      call error(this_sub_name,this_mod_name,message(1:2))
    endif
    call read_octave_al(e,'sub_e',fu) ! allocate( e )
    call read_octave_al(t,'sub_t',fu) ! allocate( t )
   close(fu,iostat=ierr)

   ! 1) count the subgrid vertices
   allocate( v2v(grid%nv) ); v2v = -1
   nv = 0 ! counter of the subgrid vertices
   do ie=1,size(t,2)
     do iv=1,grid%d ! vertices of the subgrid element
       if(v2v(t(iv,ie)).eq.-1) then ! new vertex
         nv = nv+1
         v2v(t(iv,ie)) = nv
       endif
     enddo
   enddo

   ! 2) set sub_v2v and sub_e2s (and reorder t)
   allocate( sub_v2v(nv) )
   do iv=1,grid%nv
     if(v2v(iv).ne.-1) sub_v2v(v2v(iv)) = iv
   enddo
   allocate( sub_e2s(size(t,2)) )
   do ie=1,size(t,2)
     ! we need to find the side in grid
     search_do: do isl=1,grid%v(t(1,ie))%ns
       found = .true.
       is = grid%v(t(1,ie))%is(isl)
       check_do: do iv=2,grid%d ! sides of the other vertices
         if(all(grid%v(t(iv,ie))%is.ne.is)) then
           found = .false.
           exit check_do
         endif
       enddo check_do
       if(found) exit search_do
     enddo search_do
     if(found) then
       sub_e2s(ie) = is
       t(1:grid%d,ie) = grid%s(is)%iv
     else
       write(message(1),'(a,i7,a)') &
         'Problems locating subelement ',ie,' in the original grid.'
       call error(this_sub_name,this_mod_name,message(1))
     endif
   enddo

   ! 3) set p, e, t
   allocate(p(grid%m,nv))
   do iv=1,nv
     p(:,iv) = grid%v(sub_v2v(iv))%x
   enddo
   do ib=1,size(e,2)
     do iv=1,grid%d-1
       e(iv,ib) = v2v(e(iv,ib))
     enddo
   enddo
   do ie=1,size(t,2)
     do iv=1,grid%d ! vertices of the subgrid element
       t(iv,ie) = v2v(t(iv,ie))
     enddo
   enddo

   ! 4) define the subgrid
   call new_grid( sub_grid, grid%d-1, p, e, t )

   deallocate(p,e,t,v2v)
 end subroutine new_subgrid

!-----------------------------------------------------------------------
 
 subroutine write_grid_struct(grid,var_name,fu)
 ! octave output for the complete grid
  integer, intent(in) :: fu
  type(t_grid), intent(in) :: grid
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_grid_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 12' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%d,'<cell-element>',fu)

   ! field 02 : m
   write(fu,'(a)')      '# name: m'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%m,'<cell-element>',fu)

   ! field 03 : ne
   write(fu,'(a)')      '# name: ne'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ne,'<cell-element>',fu)

   ! field 04 : nv
   write(fu,'(a)')      '# name: nv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%nv,'<cell-element>',fu)

   ! field 05 : ns
   write(fu,'(a)')      '# name: ns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ns,'<cell-element>',fu)

   ! field 06 : nb
   write(fu,'(a)')      '# name: nb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%nb,'<cell-element>',fu)

   ! field 07 : ni
   write(fu,'(a)')      '# name: ni'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ni,'<cell-element>',fu)

   ! field 08 : me
   write(fu,'(a)')      '# name: me'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%me,'<cell-element>',fu)

   ! field 09 : v
   write(fu,'(a)')      '# name: v'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%v,'<cell-element>',fu)

   ! field 10 : s
   write(fu,'(a)')      '# name: s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%s,'<cell-element>',fu)

   ! field 11 : e
   write(fu,'(a)')      '# name: e'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%e,'<cell-element>',fu)

   ! field 12 : vol
   write(fu,'(a)')      '# name: vol'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%vol,'<cell-element>',fu)

 end subroutine write_grid_struct
 
!-----------------------------------------------------------------------
 
 subroutine write_vert_struct(v,var_name,fu)
  integer, intent(in) :: fu
  type(t_v), intent(in) :: v(:)
  character(len=*), intent(in) :: var_name
 
  integer :: iv
  character(len=*), parameter :: &
    this_sub_name = 'write_vert_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 7' ! number of fields

   ! field 01 : m
   write(fu,'(a)')      '# name: m'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%m,'<cell-element>',fu)
   enddo

   ! field 02 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%i,'<cell-element>',fu)
   enddo

   ! field 03 : x
   write(fu,'(a)')      '# name: x'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%x,'c','<cell-element>',fu)
   enddo

   ! field 04 : ns
   write(fu,'(a)')      '# name: ns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%ns,'<cell-element>',fu)
   enddo

   ! field 05 : is
   write(fu,'(a)')      '# name: is'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%is,'r','<cell-element>',fu)
   enddo

   ! field 06 : ne
   write(fu,'(a)')      '# name: ne'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%ne,'<cell-element>',fu)
   enddo

   ! field 07 : ie
   write(fu,'(a)')      '# name: ie'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%ie,'r','<cell-element>',fu)
   enddo

 end subroutine write_vert_struct
 
!-----------------------------------------------------------------------
 
 subroutine write_side_struct(s,var_name,fu)
  integer, intent(in) :: fu
  type(t_s), intent(in) :: s(:)
  character(len=*), intent(in) :: var_name
 
  integer :: is
  character(len=*), parameter :: &
    this_sub_name = 'write_side_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 10' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%d,'<cell-element>',fu)
   enddo

   ! field 02 : m
   write(fu,'(a)')      '# name: m'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%m,'<cell-element>',fu)
   enddo

   ! field 03 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%i,'<cell-element>',fu)
   enddo

   ! field 04 : iv
   write(fu,'(a)')      '# name: iv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%iv,'r','<cell-element>',fu)
   enddo

   ! field 05 : ie
   write(fu,'(a)')      '# name: ie'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%ie,'r','<cell-element>',fu)
   enddo

   ! field 06 : isl
   write(fu,'(a)')      '# name: isl'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%isl,'r','<cell-element>',fu)
   enddo

   ! field 07 : xb
   write(fu,'(a)')      '# name: xb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%xb,'c','<cell-element>',fu)
   enddo

   ! field 08 : a
   write(fu,'(a)')      '# name: a'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%a,'<cell-element>',fu)
   enddo

   ! field 09 : b
   write(fu,'(a)')      '# name: b'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%b,'<cell-element>',fu)
   enddo

   ! field 10 : x0
   write(fu,'(a)')      '# name: x0'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%x0,'c','<cell-element>',fu)
   enddo

 end subroutine write_side_struct

!-----------------------------------------------------------------------

 subroutine write_elem_struct(e,var_name,fu)
  integer, intent(in) :: fu
  type(t_e), intent(in) :: e(:)
  character(len=*), intent(in) :: var_name
 
  integer :: ie
  character(len=*), parameter :: &
    this_sub_name = 'write_elem_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 16' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%d,'<cell-element>',fu)
   enddo

   ! field 02 : m
   write(fu,'(a)')      '# name: m'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%m,'<cell-element>',fu)
   enddo

   ! field 03 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%i,'<cell-element>',fu)
   enddo

   ! field 04 : iv
   write(fu,'(a)')      '# name: iv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%iv,'r','<cell-element>',fu)
   enddo

   ! field 05 : is
   write(fu,'(a)')      '# name: is'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%is,'r','<cell-element>',fu)
   enddo

   ! field 06 : pi
   write(fu,'(a)')      '# name: pi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%pi,'r','<cell-element>' , fu )
   enddo

   ! field 07 : ip
   write(fu,'(a)')      '# name: ip'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%ip,'r','<cell-element>' , fu )
   enddo

   ! field 08 : ie
   write(fu,'(a)')      '# name: ie'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%ie,'r','<cell-element>',fu)
   enddo

   ! field 09 : iel
   write(fu,'(a)')      '# name: iel'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%iel,'r','<cell-element>',fu)
   enddo

   ! field 10 : xb
   write(fu,'(a)')      '# name: xb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%xb,'c','<cell-element>',fu)
   enddo

   ! field 11 : vol
   write(fu,'(a)')      '# name: vol'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%vol,'<cell-element>',fu)
   enddo

   ! field 12 : b
   write(fu,'(a)')      '# name: b'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%b,'<cell-element>',fu)
   enddo

   ! field 13 : det_b
   write(fu,'(a)')      '# name: det_b'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%det_b,'<cell-element>',fu)
   enddo

   ! field 14 : bi
   write(fu,'(a)')      '# name: bi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%bi,'<cell-element>',fu)
   enddo

   ! field 15 : x0
   write(fu,'(a)')      '# name: x0'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%x0,'c','<cell-element>',fu)
   enddo

   ! field 16 : n
   write(fu,'(a)')      '# name: n'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%n,'<cell-element>',fu)
   enddo

 end subroutine write_elem_struct

!-----------------------------------------------------------------------

 subroutine read_al_ddc_v_struct(s,var_name,fu)
  integer, intent(in) :: fu
  type(t_ddc_nv), allocatable, intent(out) :: s(:)
  character(len=*), intent(in) :: var_name

  integer :: ierr, inrow, incol, iv
  character(len=1000) :: message
  character(len=*), parameter :: &
    this_sub_name = 'read_al_ddc_v_struct'

   call locate_var(fu,var_name,ierr)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else

     ! Read the structure definition
     read(fu,'(A)') message ! exclude "# type: struct"
     read(fu,'(A)') message ! exclude "# length: 3"

     ! Read the first field: nd
     call locate_var(fu,'nd',ierr, norewind=.true.) ! exclude "# name: nd"
     read(fu,'(A)') message ! exclude "# type: cell"
     read(fu, '(A7,I10)') message, inrow ! "# rows:"
     read(fu,'(A10,I10)') message, incol ! "# columns:"
     allocate( s(inrow*incol) )
     do iv=1,size(s)
       s(iv)%i = iv
       call read_octave(s(iv)%nd,'<cell-element>',fu, norewind=.true. )
     enddo

     ! Read the second field: id
     call locate_var(fu,'id',ierr, norewind=.true. )
     ! no need to read inrow and incol, since they are always the same
     do iv=1,size(s)
       call read_octave_al(s(iv)%id,'<cell-element>',fu, norewind=.true. )
     enddo

     ! Read the third field: in
     call locate_var(fu,'iv',ierr, norewind=.true. )
     ! no need to read inrow and incol, since they are always the same
     do iv=1,size(s)
       call read_octave_al(s(iv)%in,'<cell-element>',fu, norewind=.true. )
     enddo

   endif
  
 end subroutine read_al_ddc_v_struct

!-----------------------------------------------------------------------
 
 subroutine write_ddc_grid_struct(ddc_grid,var_name,fu)
 ! octave output for the complete grid
  integer, intent(in) :: fu
  type(t_ddc_grid), intent(in) :: ddc_grid
  character(len=*), intent(in) :: var_name
 
  integer, allocatable :: tmp(:,:)
  character(len=*), parameter :: &
    this_sub_name = 'write_ddc_grid_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 15' ! number of fields

   ! field 01 : nd
   write(fu,'(a)')      '# name: nd'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nd,'<cell-element>',fu)

   ! field 02 : id
   write(fu,'(a)')      '# name: id'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%id,'<cell-element>',fu)

   ! field 03 : ngv
   write(fu,'(a)')      '# name: ngv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%ngv,'<cell-element>',fu)

   ! field 04 : ngs
   write(fu,'(a)')      '# name: ngs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%ngs,'<cell-element>',fu)

   ! field 05 : nge
   write(fu,'(a)')      '# name: nge'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nge,'<cell-element>',fu)

   ! field 06 : nnv
   write(fu,'(a)')      '# name: nnv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nnv,'<cell-element>',fu)

   ! field 07 : nns
   write(fu,'(a)')      '# name: nns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nns,'<cell-element>',fu)

   ! field 08 : nnv_id
   write(fu,'(a)')      '# name: nnv_id'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nnv_id,'r','<cell-element>',fu)

   ! field 09 : nns_id
   write(fu,'(a)')      '# name: nns_id'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nns_id,'r','<cell-element>',fu)

   ! field 10 : ddc_marker
   write(fu,'(a)')      '# name: ddc_marker'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%ddc_marker,'<cell-element>',fu)

   ! field 11 : gv
   write(fu,'(a)')      '# name: gv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   allocate(tmp(2,size(ddc_grid%gv)))
   tmp(1,:) = ddc_grid%gv%gi
   tmp(2,:) = ddc_grid%gv%ni
   call write_octave(tmp,'<cell-element>',fu)
   deallocate(tmp)

   ! field 12 : gs
   write(fu,'(a)')      '# name: gs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   allocate(tmp(2,size(ddc_grid%gs)))
   tmp(1,:) = ddc_grid%gs%gi
   tmp(2,:) = ddc_grid%gs%ni
   call write_octave(tmp,'<cell-element>',fu)
   deallocate(tmp)

   ! field 13 : ge
   write(fu,'(a)')      '# name: ge'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%ge%gi,'r','<cell-element>',fu)

   ! field 14 : nv
   write(fu,'(a)')      '# name: nv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(ddc_grid%nv,'<cell-element>',fu)

   ! field 15 : ns
   write(fu,'(a)')      '# name: ns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   allocate(tmp(4,size(ddc_grid%ns)))
   tmp(1,:) = ddc_grid%ns%i
   tmp(2,:) = ddc_grid%ns%id
   tmp(3,:) = ddc_grid%ns%in
   tmp(4,:) = ddc_grid%ns%p_s2s
   call write_octave(tmp,'<cell-element>',fu)
   deallocate(tmp)

 end subroutine write_ddc_grid_struct
 
!-----------------------------------------------------------------------
 
 subroutine write_ddc_vert_struct(v,var_name,fu)
  integer, intent(in) :: fu
  type(t_ddc_nv), intent(in) :: v(:)
  character(len=*), intent(in) :: var_name
 
  integer :: iv
  character(len=*), parameter :: &
    this_sub_name = 'write_ddc_vert_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 4' ! number of fields

   ! field 01 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%i,'<cell-element>',fu)
   enddo

   ! field 02 : nd
   write(fu,'(a)')      '# name: nd'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%nd,'<cell-element>',fu)
   enddo

   ! field 03 : id
   write(fu,'(a)')      '# name: id'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%id,'r','<cell-element>',fu)
   enddo

   ! field 04 : in
   write(fu,'(a)')      '# name: in'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(v)
   do iv=1,size(v)
     call write_octave(v(iv)%in,'r','<cell-element>',fu)
   enddo

 end subroutine write_ddc_vert_struct
 
!-----------------------------------------------------------------------

end module mod_grid

